Matrix

The Matrix object

A Matrix in JDemetra is stored in a single array of doubles. The data are stacked column by column, following the Fortran convention. So, a 3x2 matrix (M) is internally defined as:

nrows_=3
ncols=2
data_={a00,a10,a20,a01,a11,a21}

For retrieving an individual cell, we use:

double a10=M.get(1,0);

Accessing a matrix cell by cell is expensive. The library is designed for retrieving range of cells, by means of DataBlocks.

Apart from DataBlocks, matrices can be efficiently handled by means of SubMatrix objects, which are simply partial views of the underlying matrix. It is important to note that a Matrix owns its data, while a DataBlock or a SubMatrix refers to data contained in a given Matrix: changing a DataBlock or a SubMatrix modifies the underlying Matrix.

DataBlock row0=M.row(0);
DataBlock col1=M.column(1);
DataBlock diag=M.diagonal();

Operations on matrices

The following tables contain a summary of the main methods of the (Sub)Matrix class and of other related classes. They are classified by accessors (means to retrieve parts of the data), modifiers (means to modify the state of the Matrix) and by operators (functions that generate new matrices). It should be noted that many methods are defined either for convenience or for efficiency in the Matrixclass and in the SubMatrixclass. Some methods are only defined in SubMatrix. In that case, users can transform any Matrix X in a full SubMatrix by means of the "X.subMatrix()" method. If need be, a SubMatrix S can also be transformed in a new Matrix X through the constructor "X=new Matrix(S)". However, in that case, X and S will refer to two distinct sets of data. The design only considers one type of matrix. Specialized algorithms are grouped as static methods in independent classes. Finally, it is clear that the set of methods is redundant. However, that redundancy is explained by a concern for efficiency.

Accessors

Method	Description
column(i)	X(.,i)
row(i)	X(i,.)
diagonal()	X(k,k)
subDiagonal(i)	X(k+i,k)
columns()	{X(.,i)}
rows()	{X(i,.)}
get(i,j)	X(i,j)
subMatrix (r0,r1,c0,c1)	X(i,j), r0≤i<r1 c0≤j<c1

Modifiers

Method	Description
set(i,j, x)	X(i,j)=x
add(i,j, x)	X(i,j)=x+X(i,j)
mul(i,j, x)	X(i,j)=x.X(i,j)
set(x)	∀i,j: X(i,j)=x
copy(M)	X=M
setAY(a, Y)	X=a*Y
addAY(a, Y)	X=X+ a*Y
add(M)	X=X+M
add(d)	X=X+d
sub(M)	X=X-M
sub(d)	X=X-d
mul(d)	X=X.d
chs()	X=-X
toLower()	∀i,j,i<j: X(i,j)=0
toUpper()	∀i,j,i>j: X(i,j)=0

Operators

Method	Description
plus(M)	Y=X+M
plus(d)	Y=X+d
minus(M)	Y=X-M
minus(d)	Y=X-d
times(M)	Y=X.M
times(d)	Y=X.d
transpose()	Y=X'
sum()	x=∑X(i,j)
nrm2()	x=‖X(i,j) ‖^2

Specialized algorithms

In the following tables, S, L, U, X, v indicate respectively a symmetric matrix, a lower triangular matrix, an upper triangular matrix, a general matrix and a vector (or row matrix).

SymmetricMatrix

Method	Description
fromLower(X)	∀i,j,i<j: X(i,j)=X(j,i)
fromUpper(X)	∀i,j,i>j: X(i,j)=X(j,i)
XtX(X)	Y=X'X
XXt(X)	Y=XX'
UtU(U)	Y=U'U
LLt(L)	Y=LL'
quadraticForm(S,X)	Y=X'SX
quadraticFormT(S,X)	Y=XSX'
quadraticForm(S,v)	y=vSv'
lCholesky(S)	L:LL'=S
uCholesky(S)	U:U'U=S
inverse(S)	Y=S^(-1)

LowerTriangularMatrix

Method	Description
inverse(L)	Y=L^-1
lsolve(L, X)	Y:YL=X
rsolve(L, X)	Y:LY=X
lmul(L, X)	Y=XL
rmul(L, X)	Y=LX

UpperTriangularMatrix

Method	Description
inverse(U)	Y=U^-1
lsolve(U, X)	Y:YU=X
rsolve(U, X)	Y:UY=X
lmul(U, X)	Y=XU
rmul(U, X)	Y=UX

Other advanced algorithms

ToeplitzMatrix (T)

Method	Description
inverse(T)	Y=T^(-1)
mul(T, v)	y=Tv (product by fast Fourier transform)
solveDurbinSystem
solveLevinsonSystem

Householder, HouseholderC, HouseholderR, ThinHouseholder

QR decompositions and related algorithms.

SparseSystemSolver

Resolution of a set of linear equations by means of a gaussian elimination with pivoting. Designed for sparse matrices.

Gauss, CroutDoolittle

Resolution of a set of linear equations by means of a LU decomposition with pivoting.

HouseholderReflection, GivensRotation, ElementaryTransformations

Elementary transformations of vectors.

EigenSystem

Computation of eigen values/vectors by means of the Hessenberg's algorithm.

Examples

Basic

// Matrix creation
// 5 x 5 identity matrix
Matrix I=Matrix.identity(5);
// 5 x 5 diagonal matrix
Matrix D=Matrix.diagonal(new double[]{1,2,3,4,5});
       
// 5 x 10(random) matrix
Matrix M=new Matrix (5, 10);
M.randomize();
    
// Matrix concatenations. No explicit functions
// A = I | D
Matrix A=new Matrix(5, 10);
A.subMatrix(0, 5, 0, 5).copy(I.subMatrix());
A.subMatrix(0, 5, 5, 10).copy(D.subMatrix());
    
//     M
// B = -
//     A
Matrix B=new Matrix(10, 10);
B.subMatrix(0, 5, 0, 10).copy(M.subMatrix());
B.subMatrix(5, 10, 0, 10).copy(A.subMatrix());
   
// or equivalently
// B.subMatrix(5, 10, 0, 5).copy(I.subMatrix());
// B.subMatrix(5, 10, 5, 10).copy(D.subMatrix());
    
    
// D kronecker M (not the most efficient solution)
Matrix Q=new Matrix(25, 50);
Q.subMatrix().kronecker(D.subMatrix(), M.subMatrix());
    
// Simple matrix computation
// D * M
Matrix DM = D.times(M);
// A + M
Matrix ApM = A.plus(M);
// DM * AM'
Matrix X1=DM.times(ApM.transpose());
// or better (faster)
Matrix X2=new Matrix(5, 5);
X2.subMatrix().product(DM.subMatrix(), ApM.subMatrix().transpose());
    
// Using the methods on data blocks, we can easily implements many functions
// average by columns
double[] m=new double[X1.getColumnsCount()];
    
for (int c=0; c<X1.getColumnsCount(); ++c){
   DataBlock col=X1.column(c);
   m[c]=col.sum()/col.getLength();
}

// or equivalently (if we want more special information)
double[] e=new double[X1.getColumnsCount()];
double[] med=new double[X1.getColumnsCount()];
DataBlockIterator columns=X1.columns();
DataBlock column=columns.getCurrent();
do{
   DescriptiveStatistics stats=new DescriptiveStatistics(column);
   m[columns.getPosition()]=stats.getAverage();
   e[columns.getPosition()]=stats.getStdev();
   med[columns.getPosition()]=stats.getMedian();
}while (columns.next());

Advanced

The following code shows how JTsToolkit can be used to solve the common problem:

Ax = B

It should be noted that users should nearly never solve that problem by computing

x=A^(-1) B

Creates a 100 x 100 random matrix (A) and a 100 x 20 random matrix (B) 
int n = 100, m=20;
Matrix A = new Matrix(n, n);
A.randomize();
Matrix B = new Matrix(n, 20);
B.randomize();

// Use Householder (without pivoting). A good compromise between speed and accuracy.
ILinearSystemSolver solver = new Householder(false);
solver.decompose(A.clone());
Matrix C1 = solver.solve(B);

// Use LU decomposition (with pivoting). The fastest solution of the library (20% faster than Householder).
ILinearSystemSolver solver = new CroutDoolittle();
solver.decompose(A.clone());
Matrix C2 = solver.solve(B);