Matrix Decompositions in IronPython QuickStart Sample

Illustrates how compute various decompositions of a matrix using classes in the Numerics.NET.LinearAlgebra namespace in IronPython.

View this sample in: C# Visual Basic F#

```Python
import numerics

from System import Array

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

# Illustrates the use of matrix decompositions for solving systems of
# simultaneous linear equations and related operations using the 
# Decomposition class and its derived classes from the
# Extreme.Mathematics.LinearAlgebra namespace of Numerics.NET.

# For details on the basic workings of Vector 
# objects, including constructing, copying and
# cloning vectors, see the BasicVectors QuickStart
# Sample.
#
# For details on the basic workings of Matrix
# objects, including constructing, copying and
# cloning vectors, see the BasicVectors QuickStart
# Sample.
#

#
# LU Decomposition
#

# The LU decomposition of a matrix rewrites a matrix A in the
# form A = PLU with P a permutation matrix, L a unit-
# lower triangular matrix, and U an upper triangular matrix.

aLU = Matrix([ \
    [ 1.80, 2.88, 2.05,-0.89 ], \
    [ 5.25,-2.95,-0.95,-3.80 ], \
    [ 1.58,-2.69,-2.90,-1.04 ], \
    [-1.11,-0.66,-0.59, 0.80 ]])

bLU = Matrix([[ 9.52, 18.47], [24.35,2.25,], [0.77,-13.28], [-6.22,-6.21 ]])

# The decomposition is obtained by calling the GetLUDecomposition
# method of the matrix. It takes zero or one parameters. The
# parameter is a bool value that indicates whether the
# matrix may be overwritten with its decomposition.
lu = aLU.GetLUDecomposition(True)
print "A = {0:F2}", aLU

# The Decompose method performs the decomposition. You don't need
# to call it explicitly, as it is called automatically as needed.
			
# The IsSingular method checks for singularity.
print "'A is singular' is", lu.IsSingular()
# The LowerTriangularFactor and UpperTriangularFactor properties
# return the two main components of the decomposition.
print "L = {0:.6f}.".format(lu.LowerTriangularFactor)
print "U = {0:.6f}.".format(lu.UpperTriangularFactor)

# GetInverse() gives the matrix inverse, Determinant() the determinant:
print "Inv A = {0:.6f}.".format(lu.GetInverse())
print "Det A = {0:.6f}.".format(lu.GetDeterminant())

# The Solve method solves a system of simultaneous linear equations, with
# one or more right-hand-sides:
xLU = lu.Solve(bLU)
print "x = {0:.6f}.".format(xLU)

# The permutation is available through the RowPermutation property:
print "P =",  lu.RowPermutation
print "Px = {0:.6f}.".format(xLU.PermuteRows(lu.RowPermutation))

#
# QR Decomposition
#

# The QR decomposition of a matrix A rewrites the matrix
# in the form A = QR, with Q a square, orthogonal matrix, # and R an upper triangular matrix.

aQR = Matrix([[2.0,2.5,2.5], [2.0,2.5,2.5], [1.6,-0.4,2.8], [2.0,-0.5,0.5],[1.2,-0.3,-2.9]])
bQR = Vector([1.1, 0.9, 0.6, 0.0,-0.8])

# The decomposition is obtained by calling the GetQRDecomposition
# method of the matrix. It takes zero or one parameters. The
# parameter is a bool value that indicates whether the
# matrix may be overwritten with its decomposition.
qr = aQR.GetQRDecomposition(True)
print "A = {0:.1f}.".format(aQR)

# The Decompose method performs the decomposition. You don't need
# to call it explicitly, as it is called automatically as needed.
			
# The IsSingular method checks for singularity.
print "'A is singular' is", qr.IsSingular()

# GetInverse() gives the matrix inverse, Determinant() the determinant, # but these are defined only for square matrices.

# The Solve method solves a system of simultaneous linear equations, with
# one or more right-hand-sides. If the matrix is over-determined, you can
# use the LeastSquaresSolve method to find a least squares solution:
xQR = qr.LeastSquaresSolve(bQR)
print "x = {0:.6f}.".format(xQR)

# The OrthogonalFactor and UpperTriangularFactor properties
# return the two main components of the decomposition.
print "Q = {0:.6f}.".format(qr.OrthogonalFactor.AsDenseMatrix())
print "R = {0:.6f}.".format(qr.UpperTriangularFactor)

# You don't usually need to form Q explicitly. You can multiply
# a vector or matrix on either side by Q using the Multiply method:
print "Qx = {0:.6f}.".format(qr.OrthogonalFactor.Multiply(xQR))
print "transpose(Q)x = {0:.6f}.".format(qr.OrthogonalFactor.MultiplyTranspose(xQR))

#
# Singular Value Decomposition
#

# The singular value decomposition of a matrix A rewrites the matrix
# in the form A = USVt, with U and V orthogonal matrices, # S a diagonal matrix. The diagonal elements of S are called 
# the singular values.

aSvd = Matrix([[2.0,2.0,1.6,2.0,1.2],[2.5,2.5,-0.4,-0.5,-0.3],[2.5, 2.5, 2.8, 0.5,-2.9]])
bSvd = Vector([1.1, 0.9, 0.6])

# The decomposition is obtained by calling the GetSingularValueDecomposition
# method of the matrix. It takes zero or one parameters. The
# parameter indicates which parts of the decomposition
# are to be calculated. The default is All.
svd = aSvd.GetSingularValueDecomposition()
print "A = {0:.1f}.".format(aSvd)

# The Decompose method performs the decomposition. You don't need
# to call it explicitly, as it is called automatically as needed.
			
# The IsSingular method checks for singularity.
print "'A is singular' is", svd.IsSingular()

# Several methods are specific to this class. The GetPseudoInverse
# method returns the Moore-Penrose pseudo-inverse, a generalization
# of the inverse of a matrix to rectangular and/or singular matrices:
aInv = svd.GetPseudoInverse()

# It can be used to solve over- or under-determined systems.
xSvd = aInv * bSvd
print "x = {0:.6f}.".format(xSvd)

# The SingularValues property returns a vector that contains
# the singular values in descending order:
print "S = {0:.6f}".format(svd.SingularValues)

# The LeftSingularVectors and RightSingularVectors properties
# return matrices that contain the U and V factors
# of the decomposition.
print "U = {0:.6f}.".format(svd.LeftSingularVectors)
print "V = {0:.6f}.".format(svd.RightSingularVectors)

#
# Cholesky decomposition
#

# The Cholesky decomposition of a symmetric matrix A
# rewrites the matrix in the form A = GGt with
# G a lower-triangular matrix.

# Remember the column-major storage mode: each line of
# components contains one COLUMN of the matrix.
aC = Matrix.CreateSymmetric(4, Array[float]([ \
	4.16,-3.12, 0.56,-0.10, \
    0, 5.03,-0.83, 1.18, \
    0,0, 0.76, 0.34, \
    0,0,0, 1.18 ]), \
    MatrixTriangle.Lower)
bC = Matrix([[8.70,8.30], [-13.35,2.13], [1.89,1.61], [-4.14,5.00]])

# The decomposition is obtained by calling the GetCholeskyDecomposition
# method of the matrix. It takes zero or one parameters. The
# parameter is a bool value that indicates whether the
# matrix should be overwritten with its decomposition.
c = aC.GetCholeskyDecomposition(True)
print "A = {0:F2}", aC

# The Decompose method performs the decomposition. You don't need
# to call it explicitly, as it is called automatically as needed.
			
# The IsSingular method checks for singularity.
print "'A is singular' is", c.IsSingular()
# The LowerTriangularFactor returns the component of the decomposition.
print "L = {0:.6f}.".format(c.LowerTriangularFactor)

# GetInverse() gives the matrix inverse, Determinant() the determinant:
print "Inv A = {0:.6f}.".format(c.GetInverse())
print "Det A = {0:.6f}.".format(c.GetDeterminant())

# The Solve method solves a system of simultaneous linear equations, with
# one or more right-hand-sides:
xC = c.Solve(bC)
print "x = {0:.6f}.".format(xC)
			
#
# Symmetric eigenvalue decomposition
#

# The eigenvalue decomposition of a symmetric matrix A
# rewrites the matrix in the form A = XLXt with
# X an orthogonal matrix and L a diagonal matrix.
# The diagonal elements of L are the eigenvalues.
# The columns of X are the eigenvectors.

# Remember the column-major storage mode: each line of
# components contains one COLUMN of the matrix.
aEig = Matrix.CreateSymmetric(4, Array[float]([ \
	0.5,  0.0,  2.3, -2.6, \
    0.0,  0.5, -1.4, -0.7, \
    2.3, -1.4,  0.5,  0.0, \
   -2.6, -0.7,  0.0,  0.5]), MatrixTriangle.Lower)

# The decomposition is obtained by calling the GetLUDecomposition
# method of the matrix. It takes zero or one parameters. The
# parameter is a bool value that indicates whether the
# matrix should be overwritten with its decomposition.
eig = aEig.GetEigenvalueDecomposition()
print "A = {0:.2f}.".format(aEig)

# The Decompose method performs the decomposition. You don't need
# to call it explicitly, as it is called automatically as needed.
			
# The IsSingular method checks for singularity.
print "'A is singular' is", eig.IsSingular()
# The eigenvalues and eigenvectors of a symmetric matrix are all real.
# The RealEigenvalues property returns a vector containing the eigenvalues:
print "L = {0:.6f}.".format(eig.RealEigenvalues)
# The RealEigenvectors property returns a matrix whose columns
# contain the corresponding eigenvectors:
print "X = {0:.6f}.".format(eig.RealEigenvectors)

# GetInverse() gives the matrix inverse, Determinant() the determinant:
print "Inv A = {0:.6f}.".format(eig.GetInverse())
print "Det A = {0:.6f}.".format(eig.GetDeterminant())
```