Illustrates how compute various decompositions of a matrix using classes in the Extreme.Mathematics.LinearAlgebra namespace in IronPython.
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 the Extreme Optimization
#/ Mathematics Library for .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())