Least Squares in IronPython QuickStart Sample

Illustrates how to solve least squares problems using classes in the Extreme.Mathematics.LinearAlgebra namespace in IronPython.

View this sample in: C# Visual Basic F#

```Python
import numerics

# The DenseMatrix and LeastSquaresSolver classes reside in the 
# Extreme.Mathematics.LinearAlgebra namespace.
from Extreme.Mathematics import *
from Extreme.Mathematics.LinearAlgebra import *

# Illustrates solving least squares problems using the 
# LeastSquaresSolver class in the Extreme.Mathematics.LinearAlgebra 
# namespace of Extreme Numerics.NET.

# A least squares problem consists in finding
# the solution to an overdetermined system of
# simultaneous linear equations so that the
# sum of the squares of the error is minimal.
#
# A common application is fitting data to a
# curve. See the CurveFitting sample application
# for a complete example.

# Let's start with a general matrix. This will be
# the matrix a in the left hand side ax=b:
a = Matrix([[1,1,1,1], [1,2,4,2],[1,3,9,1],[1,4,16,2],[1,5,25,1],[1,6,36,2]])
# Here is the right hand side:
b = Vector([1, 3, 6, 11, 15, 21])
b2 = Matrix([[1,1],[3,2],[6,3],[11,4],[15,5],[21,7]])
print "a = {0:.4f}".format(a)
print "b = {0:.4f}".format(b)

#
# The LeastSquaresSolver class
#

# The following creates an instance of the
# LeastSquaresSolver class for our problem:
solver = LeastSquaresSolver(a, b)
# We can specify the solution method: normal
# equations or QR decomposition. In most cases, # a QR decomposition is the most desirable:
solver.SolutionMethod = LeastSquaresSolutionMethod.QRDecomposition
# The Solve method calculates the solution:
x = solver.Solve()
print "x = {0:.4f}".format(x)
# The Solution property also returns the solution:
print "x = {0:.4f}".format(solver.Solution)
# More detailed information is available from
# additional methods.
# The values of the right hand side predicted 
# by the solution:
print "Predictions = {0:.4f}".format(solver.GetPredictions())
# The residuals (errors) of the solution:
print "Residuals = {0:.4f}".format(solver.GetResiduals())
# The total sum of squares of the residues:
print "Residual square error =", solver.GetResidualSumOfSquares()

#
# Direct normal equations
#

# Alternatively, you can create a least squares
# solution by providing the normal equations
# directly. This may be useful when it is easy
# to calculate the normal equations directly.
# 
# Here, we'll just calculate the normal equation:
aTa = SymmetricMatrix.FromOuterProduct(a)
aTb = b * a # a.Transpose() * b
# We find the solution by solving the normal equations
# directly:
x = aTa.Solve(aTb)
print "x = {0:.4f}".format(x)
# However, properties of the least squares solution, such as
# error estimates and residuals are not available.
```