New Version 6.0!

Try it for free with our fully functional 60-day trial version.

Download now!

QuickStart Samples

Matrix Decompositions QuickStart Sample (C#)

Illustrates how compute various decompositions of a matrix using classes in the Extreme.Mathematics.LinearAlgebra namespace in C#.

Visual Basic code F# code IronPython code Back to QuickStart Samples

using System;

namespace Extreme.Numerics.QuickStart.CSharp
{
    using Extreme.Mathematics;

    /// <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 the Extreme Optimization
    /// Mathematics Library for .NET.
    /// </summary>
    class MatrixDecompositions
    {
        /// <summary>
        /// The main entry point for the application.
        /// </summary>
        [STAThread]
        static void Main(string[] args)
        {
            // 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();
        }
    }
}