20.3 Using Factorizations (.NET, C#, CSharp, VB, Visual Basic, F#)
Once a factorization is constructed from a matrix (see Section 20.2), it can be reused to solve for different right hand sides, and to compute inverses, determinants, and condition numbers. Static methods are also provided on MatrixFunctions to perform these functions without having to explicitly construct a factorization object.
You can use a factorization to solve for right-hand sides using the Solve() method. For instance, this code solves for one right-hand side:
Code Example – C# matrix factorization
var genMat = new DoubleMatrix(
"5x5 [ 1.0000 0.5000 0.2500 0.1250 0.0625
0.5000 1.0000 0.5000 0.2500 0.1250
0.2500 0.5000 1.0000 0.5000 0.2500
0.1250 0.2500 0.5000 1.0000 0.5000
0.0625 0.1250 0.2500 0.5000 1.0000 ]" );
var A = new DoubleSymmetricMatrix( genMat );
var F = new DoubleSymPDFact( A );
var v = new DoubleVector( A.Order, new RandGenUniform(-1,1) );
DoubleVector x = F.Solve( v );
Code Example – VB matrix factorization
Dim GenMat As New DoubleMatrix(
"5x5 [ 1.0000 0.5000 0.2500 0.1250 0.0625" &
"0.5000 1.0000 0.5000 0.2500 0.1250" &
"0.2500 0.5000 1.0000 0.5000 0.2500" &
"0.1250 0.2500 0.5000 1.0000 0.5000" &
"0.0625 0.1250 0.2500 0.5000 1.0000 ]")
Dim A As New DoubleSymmetricMatrix(GenMat)
Dim F As New DoubleSymPDFact(A)
Dim V As New DoubleVector(A.Order, New RandGenUniform(-1.0, 1.0))
Dim X As DoubleVector = F.Solve(V)
The returned vector x is the solution to the linear system Ax = v. Note that the length of vector v must be equal to the number of columns in the factored matrix A or a MismatchedSizeException is thrown.
To do the same thing without explicitly constructing a factorization object, you could do this:
Code Example – C# matrix factorization
DoubleVector x = MatrixFunctions.Solve( A, v, true );
Code Example – VB matrix factorization
Dim X As DoubleVector = MatrixFunctions.Solve(A, V, True)
The optional third, boolean parameter indicates that A is positive definite.
Similarly, you can use the Solve() method to solve for multiple right-hand sides. This code solves for 10 right-hand sides:
Code Example – C# matrix factorization
int rows = 8, cols = 8;
DoubleComplexVector data =
new DoubleComplexVector( cols*3, new RandGenUniform(-1, 1) );
DoubleComplexTriDiagMatrix A =
new DoubleComplexTriDiagMatrix( data, rows, cols );
var F = new DoubleComplexTriDiagFact( A );
var B =
new DoubleComplexMatrix( A.Cols, 10, new RandGenUniform(-1,1) );
DoubleComplexMatrix X = F.Solve( B );
Code Example – VB matrix factorization
Dim Rows As Integer = 8
Dim Cols As Integer = 8
Dim Data As New DoubleComplexVector(Cols * 3,
New RandGenUniform(-1.0, 1.0))
Dim A As New DoubleComplexTriDiagMatrix(Data, Rows, Cols)
Dim F As New DoubleComplexTriDiagFact(A)
Dim B As New DoubleComplexMatrix(A.Cols, 10,
New RandGenUniform(-1.0, 1.0))
Dim X As DoubleComplexMatrix = F.Solve(B)
The returned matrix X is the solution to the linear system AX= B. That is, the right-hand sides are the columns of B, and the solutions are the columns of X. Matrix B must have the same number of columns as the factored matrix A.
To do the same thing without explicitly constructing a factorization object, you could do this:
Code Example – C# matrix factorization
DoubleComplexMatrix X = MatrixFunctions.Solve( A, B );
Code Example – VB matrix factorization
Dim X As DoubleComplexMatrix = MatrixFunctions.Solve(A, B)
Computing Inverses, Determinants, and Condition Numbers
You can use a factorization to compute inverses using the Inverse() method, and determinants using the Determinant() method. For example:
Code Example – C# matrix factorization
int rows = 8, cols = 8;
var Lehmer = new FloatComplexMatrix( rows, cols );
for ( int i = 0; i < rows; ++i )
{
for ( int j = 0; j < cols; ++j )
{
if ( j >= i )
{
Lehmer[i,j] = (float)(i+1)/(float)(j+1);
}
}
}
var A = new FloatHermitianMatrix( Lehmer );
var F = new FloatHermitianPDFact( A );
FloatHermitianMatrix AInv = F.Inverse();
FloatComplex det = F.Determinant();
Code Example – VB matrix factorization
Dim Rows As Integer = 8
Dim Cols As Integer = 8
Dim Lehmer As New FloatComplexMatrix(Rows, Cols)
For I As Integer = 0 To Rows - 1
For J As Integer = 0 To Cols - 1
If J >= I Then
Lehmer(I, J) = CType(I + 1, Single) / CType(J + 1, Single)
End If
Next
Next
Dim A As New FloatHermitianMatrix(Lehmer)
Dim F As New FloatHermitianPDFact(A)
Dim AInv As FloatHermitianMatrix = F.Inverse()
Dim Det As FloatComplex = F.Determinant()
The ConditionNumber() method computes an estimate of the condition number in the one-norm. The condition number of a matrix A is:
kappa = ||A|| ||AInv||
where AInv is the inverse of the matrix A. For instance:
Code Example – C# matrix factorization
DoubleMatrix genMat =
new DoubleMatrix( 25, 25, new RandGenUniform( 0, 100 ) );
var A = new DoubleSymmetricMatrix( genMat );
var F = new DoubleSymFact( A );
double cond = F.ConditionNumber();
Code Example – VB matrix factorization
Dim GenMat As New DoubleMatrix(25, 25,
New RandGenUniform(0.0, 100.0))
Dim A As New DoubleSymmetricMatrix(GenMat)
Dim F As New DoubleSymFact(A)
Dim Cond As Double = F.ConditionNumber()
NOTE—The ConditionNumber() method returns the reciprocal of the condition number, rho, where rho = 1/kappa.
Banded and tridiagonal factorization classes also provide an overload of the ConditionNumber() method that accepts a value from the NormType enumeration for specifying the matrix norm.
Thus, this code estimates the condition number in the infinity-norm:
Code Example – C# matrix factorization
int rows = 4, cols = 4;
FloatVector data =
new FloatVector( cols*3, new RandGenUniform(-1, 1) );
var A = new FloatTriDiagMatrix( data, rows, cols );
var F = new FloatTriDiagFact( A );
float cond = F.ConditionNumber( NormType.InfinityNorm );
Code Example – VB matrix factorization
Dim Rows As Integer = 4
Dim Cols As Integer = 4
Dim Data As New FloatVector(Cols * 3,
New RandGenUniform(-1.0, 1.0))
Dim A As New FloatTriDiagMatrix(Data, Rows, Cols)
Dim F As New FloatTriDiagFact(A)
Dim Cond As Single = F.ConditionNumber(NormType.InfinityNorm)
Inverses, determinants, and condition numbers can also be computed without explicitly constructing a factorization object by using static methods on MatrixFunctions. For instance:
Code Example – C# matrix factorization
DoubleMatrix genMat =
new DoubleMatrix( 12, 12, new RandGenUniform( -1, 1 ) );
var A = new DoubleSymmetricMatrix( genMat );
DoubleSymmetricMatrix AInv = MatrixFunctions.Inverse( A );
double det = MatrixFunctions.Determinant( A );
double cond = MatrixFunctions.ConditionNumber( A );
Code Example – VB matrix factorization
Dim GenMat As New DoubleMatrix(12, 12,
New RandGenUniform(-1.0, 1.0))
Dim A As New DoubleSymmetricMatrix(GenMat)
Dim AInv As DoubleSymmetricMatrix = MatrixFunctions.Inverse(A)
Dim Det As Double = MatrixFunctions.Determinant(A)
Dim Cond As Double = MatrixFunctions.ConditionNumber(A)