4.3 Using Factorizations

Once a factorization is constructed from a matrix (see Section 4.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.

Solving for Right-Hand Sides

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:

DoubleMatrix 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 ]" );
DoubleSymmetricMatrix A = new DoubleSymmetricMatrix( genMat );
DoubleSymPDFact F = new DoubleSymPDFact( A );
DoubleVector v =
   new DoubleVector( A.Order, new RandGenUniform(-1,1) );
DoubleVector x = 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:

DoubleVector x = 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:

int rows = 8, cols = 8;
DoubleComplexVector data =
  new DoubleComplexVector( cols*3, new RandGenUniform(-1, 1) );
DoubleComplexTriDiagMatrix A =
  new DoubleComplexTriDiagMatrix( data, rows, cols );
DoubleComplexTriDiagFact F = new DoubleComplexTriDiagFact( A );

DoubleComplexMatrix B =
  new DoubleComplexMatrix( A.Cols, 10, new RandGenUniform(-1,1) );
      
DoubleComplexMatrix X = 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:

DoubleComplexMatrix X = 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:

int rows = 8, cols = 8;
FloatComplexMatrix 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);
    }
  }
}
FloatHermitianMatrix A = new FloatHermitianMatrix( Lehmer );

FloatHermitianPDFact F = new FloatHermitianPDFact( A );
FloatHermitianMatrix AInv = F.Inverse();
FloatComplex det = 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:

DoubleMatrix genMat =
  new DoubleMatrix( 25, 25, new RandGenUniform( 0, 100 );
DoubleSymmetricMatrix A = new DoubleSymmetricMatrix( genMat );

DoubleSymFact F = new DoubleSymFact( A );
double cond = 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:

int rows = 4, cols = 4;
FloatVector data =
   new FloatVector( cols*3, new RandGenUniform(-1, 1) );
FloatTriDiagMatrix A = new FloatTriDiagMatrix( data, rows, cols );

FloatTriDiagFact F = new FloatTriDiagFact( A );
float cond = 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:

DoubleMatrix genMat =
  new DoubleMatrix( 12, 12, new RandGenUniform( -1, 1 );
DoubleSymmetricMatrix A = new DoubleSymmetricMatrix( genMat );

DoubleSymmetricMatrix AInv = MatrixFunctions.Inverse( A );
double det = MatrixFunctions.Determinant( A );
double cond = MatrixFunctions.ConditionNumber( A );

TOC |  Previous |  Next |  Index