# 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()) ```