Linear Systems

Back to: Main contents.

Contents

Inverses and linear systems
The Solve suite
In many applications one must solve linear systems of the form
AX = B,
where A is a nonsingular matrix. The Solve suite contains static methods for solving this and related systems. Before describing these methods, we will first discuss the use of inverses to solve systems.

Inverses and linear systems

Back to: Main contents, Top of section.

The solution of the system AX = B can be written in the form X= A-1B. This suggests the following algorithm.

1. Compute C = A-1

2. Compute X = CB
This algorithm can be implement in Jampack as follows. 
Zmat C = Inv.o(A);
Zmat X  = Times.o(C, B);
Or more succinctly: 
Zmat X  = Times.o(Inv.o(A), B);
On the other hand, we can also solve the problem using the function aib from the Solve suite. 
Zmat X = Solve.aib(A, B);
The question is: Which is better?

Without going into details, the method used by Solve is better for two reasons.

1. It is faster.

2. It is numerically more stable.
For these reasons we recommend that the Jampack user substitute an appropriate Solve method for an explicit inverse wherever possible. For example, suppose that we have a formula of the form
S = A22 - A12*A11-1*A12.
Instead of computing an inverse and using the standard matrix functions to compute the matrix S, we can recognize that the matrix A11-1*A12 is the solution of the system A11*X=A12. Hence if we write
S = A22 - A12*X,
we have the following algorithm 
Zmat X = Solve.aib(A11, A12);
Zmat S = Minus.o(A22, Times.o(A12, X));
In most instances when a matrix inverse appears in a formula, it can be replaced with an appropriate Solve method.

The Solve suite

Back to: Main contents, Top of section.

In spite of the large number of individual routines in Solve, the class methods can be described in a couple of tables. (For a complete listing of the calling sequences see the Java documentation or the raw code.)

First, the Solve suite solves four kinds of linear systems involving a nonsingular matrix A. The following table gives the names of the Solve methods, the systems they solve, and where the names come from.

aib(A, B) AX = B A Inverse B
ahib(A, B) AHX = B AH Inverse B
bai(B, A) XA = B B Inverse A
bahi(B, A) XAH = B B AH Inverse

(Ironically, although we have deprecated the use of inverses in solving systems, we have used them to generate the names for our methods. A foolish consistency is the hobgoblin of petty minds.)

The second table lists the classes that each of four functions supports. The matrix B is always a Zmat, so we need only list the possible classes of A.

aib Zmat Zltmat Zutmat Zpsdmat
ahib Zmat Zltmat Zutmat
bai Zmat Zltmat Zutmat Zpsdmat
bahi Zmat Zltmat Zutmat

The methods ahib and bahi are not implemented Zpsdmats. This is natural, since positive definite matrices are Hermitian and the two routines are therefore redundant. Note that it is perfectly possible to call ahib and bahi with a Zpsdmat for the first argument; however, it will be treated as if the argument were a Zmat.

In all the Solve methods the user is responsible for insuring that the matrix in question is nonsingular. If the matrix is singular or nearly singular, the methods may throw a JampackException. Owing to rounding error, however, there is no guarantee that they will. A later version of Jampack will implement condition estimators to help the user decide on the quality of a system.

The Solve suite shows the powers of tag classes. Triangular systems can be solved much faster than general systems. Consequently, if a triangular matrix is properly tagged, Solve will automatically take advantage of its triangularity.

Solve also takes full advantage of the history feature. Specifically, to solve a general system, the Solve routines require a Zludpp decomposition, which is expensive to compute. However, once it has been computed, it is put on the history list and can be reused. Thus in the sequence 

while (!converged){
   x = Solve.aib(A, x);
   x = Times.o(1/Norm.fro(x), x);
   ...
}
the first iteration of the loop will be expensive but the subsequent interations will be cheap (unless A is altered in the latter part of the loop).