Matrix Decompositions in C# QuickStart Sample
Illustrates how compute various decompositions of a matrix using classes in the Extreme.Mathematics.LinearAlgebra namespace in C#.
View this sample in: Visual Basic F# IronPython
using System;
using Extreme.Mathematics;
namespace Extreme.Numerics.QuickStart.CSharp
{
/// <summary>
/// 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 Extreme Numerics.NET.
/// </summary>
class MatrixDecompositions
{
static void Main(string[] args)
{
var success = Extreme.License.Verify("Demo license");
Console.WriteLine(success);
// For details on the basic workings of var
// 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.
var aLU = Matrix.Create(4, 4,
new double[]
{
1.80, 5.25, 1.58, -1.11,
2.88,-2.95, -2.69, -0.66,
2.05,-0.95, -2.90, -0.59,
-0.89,-3.80,-1.04, 0.80
}, MatrixElementOrder.ColumnMajor);
var bLU = Matrix.Create(4, 2, new double[]
{
9.52,24.35,0.77,-6.22,
18.47,2.25,-13.28,-6.21
}, MatrixElementOrder.ColumnMajor);
// 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.
var lu = aLU.GetLUDecomposition(true);
Console.WriteLine("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.
Console.WriteLine("'A is singular' is {0:F6}.", lu.IsSingular());
// The LowerTriangularFactor and UpperTriangularFactor properties
// return the two main components of the decomposition.
Console.WriteLine("L = {0:F6}", lu.LowerTriangularFactor);
Console.WriteLine("U = {0:F6}", lu.UpperTriangularFactor);
// GetInverse() gives the matrix inverse, Determinant() the determinant:
Console.WriteLine("Inv A = {0:F6}", lu.GetInverse());
Console.WriteLine("Det A = {0:F6}", lu.GetDeterminant());
// The Solve method solves a system of simultaneous linear equations, with
// one or more right-hand-sides:
var xLU = lu.Solve(bLU);
Console.WriteLine("x = {0:F6}", xLU);
// The permutation is available through the RowPermutation property:
Console.WriteLine("P = {0}", lu.RowPermutation);
Console.WriteLine("Px = {0:F6}", xLU.PermuteRowsInPlace(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.
var aQR = Matrix.Create(5, 3,
new double[]
{
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
}, MatrixElementOrder.ColumnMajor);
var bQR = Vector.Create(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.
var qr = aQR.GetQRDecomposition(true);
Console.WriteLine("A = {0:F1}", 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.
Console.WriteLine("'A is singular' is {0:F6}.", 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:
var xQR = qr.LeastSquaresSolve(bQR);
Console.WriteLine("x = {0:F6}", xQR);
// The OrthogonalFactor and UpperTriangularFactor properties
// return the two main components of the decomposition.
Console.WriteLine("Q = {0:F6}", qr.OrthogonalFactor.ToDenseMatrix());
Console.WriteLine("R = {0:F6}", 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:
var Q = qr.OrthogonalFactor;
Console.WriteLine("Qb = {0:F6}", qr.OrthogonalFactor * bQR);
Console.WriteLine("transpose(Q)b = {0:F6}",
qr.OrthogonalFactor.Transpose() * bQR);
//
// 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.
var aSvd = Matrix.Create(3, 5,
new double[]
{
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
}, MatrixElementOrder.RowMajor);
var bSvd = Vector.Create(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.
var svd = aSvd.GetSingularValueDecomposition();
Console.WriteLine("A = {0:F1}", 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.
Console.WriteLine("'A is singular' is {0:F6}.", 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:
var aInv = svd.GetPseudoInverse();
// It can be used to solve over- or under-determined systems.
var xSvd = aInv * bSvd;
Console.WriteLine("x = {0:F6}", xSvd);
// The SingularValues property returns a vector that contains
// the singular values in descending order:
Console.WriteLine("S = {0:F6}", svd.SingularValues);
// The LeftSingularVectors and RightSingularVectors properties
// return matrices that contain the U and V factors
// of the decomposition.
Console.WriteLine("U = {0:F6}", svd.LeftSingularVectors);
Console.WriteLine("V = {0:F6}", 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.
var aC = Matrix.CreateSymmetric(4,
new double[]
{
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, MatrixElementOrder.ColumnMajor);
var bC = Matrix.Create(4, 2,
new double[] {8.70,-13.35,1.89,-4.14,8.30,2.13,1.61,5.00},
MatrixElementOrder.ColumnMajor);
// 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.
var c = aC.GetCholeskyDecomposition(true);
Console.WriteLine("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.
Console.WriteLine("'A is singular' is {0:F6}.", c.IsSingular());
// The LowerTriangularFactor returns the component of the decomposition.
Console.WriteLine("L = {0:F6}", c.LowerTriangularFactor);
// GetInverse() gives the matrix inverse, Determinant() the determinant:
Console.WriteLine("Inv A = {0:F6}", c.GetInverse());
Console.WriteLine("Det A = {0:F6}", c.GetDeterminant());
// The Solve method solves a system of simultaneous linear equations, with
// one or more right-hand-sides:
var xC = c.Solve(bC);
Console.WriteLine("x = {0:F6}", 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.
var aEig = Matrix.CreateSymmetric(4,
new double[]
{
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, MatrixElementOrder.ColumnMajor);
// 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.
var eig = aEig.GetEigenvalueDecomposition();
Console.WriteLine("A = {0:F2}", 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.
Console.WriteLine("'A is singular' is {0:F6}.", eig.IsSingular());
// The eigenvalues and eigenvectors of a symmetric matrix are all real.
// The RealEigenvalues property returns a vector containing the eigenvalues:
Console.WriteLine("L = {0:F6}", eig.Eigenvalues);
// The RealEigenvectors property returns a matrix whose columns
// contain the corresponding eigenvectors:
Console.WriteLine("X = {0:F6}", eig.Eigenvectors);
// GetInverse() gives the matrix inverse, Determinant() the determinant:
Console.WriteLine("Inv A = {0:F6}", eig.GetInverse());
Console.WriteLine("Det A = {0:F6}", eig.GetDeterminant());
Console.Write("Press Enter key to exit...");
Console.ReadLine();
}
}
}