/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.util;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.CholeskyDecomposition;
import cern.colt.matrix.linalg.SingularValueDecomposition;
/**
* This class implements a few methods for matrix calculations
* with double numbers.
*
*/
public class DMatrix {
private double [][] mat; // matrix of double's
private int r, c; // number of rows, columns
/**
* Creates a new DMatrix with r rows and
* c columns.
*
* @param r the number of rows
*
* @param c the number of columns
*
*
*/
public DMatrix (int r, int c) {
mat = new double[r][c];
this.r = r;
this.c = c;
}
/**
* Creates a new DMatrix with r rows and
* c columns using the data in data.
*
* @param data the data of the new DMatrix
*
* @param r the number of rows
*
* @param c the number of columns
*
*
*/
public DMatrix (double[][] data, int r, int c) {
this (r, c);
for(int i = 0; i < r; i++)
for(int j = 0; j < c; j++)
mat[i][j] = data[i][j];
}
/**
* Copy constructor.
*
* @param that the DMatrix to copy
*
*/
public DMatrix (DMatrix that) {
this (that.mat, that.r, that.c);
}
/**
* Given a symmetric positive-definite matrix M, performs the Cholesky
* decomposition of M and returns the result as a lower triangular matrix L,
* such that M = LLT.
*
* @param M the input matrix
*
* @param L the Cholesky lower triangular matrix
*
*
*/
public static void CholeskyDecompose (double[][] M, double[][] L) {
int d = M.length;
DoubleMatrix2D myM = new DenseDoubleMatrix2D (M);
DoubleMatrix2D myL = new DenseDoubleMatrix2D (d, d);
CholeskyDecomposition decomp = new CholeskyDecomposition (myM);
myL = decomp.getL();
for(int i = 0; i < L.length; i++)
for(int j = 0; j <= i; j++)
L[i][j] = myL.get(i,j);
for(int i = 0; i < L.length; i++)
for(int j = i + 1; j < L.length; j++)
L[i][j] = 0.0;
}
/**
* Given a symmetric positive-definite matrix M, performs the Cholesky
* decomposition of M and returns the result as a lower triangular matrix L,
* such that M = LLT.
*
* @param M the input matrix
*
* @return the Cholesky lower triangular matrix
*
*/
public static DoubleMatrix2D CholeskyDecompose (DoubleMatrix2D M) {
CholeskyDecomposition decomp = new CholeskyDecomposition (M);
return decomp.getL();
}
/**
* Computes the principal components decomposition M =
*
* UΛUt by using the singular value decomposition of matrix
* M. Puts the eigenvalues, which are the diagonal elements of matrix Λ,
* sorted by decreasing size, in vector lambda, and puts matrix
*
* A = U(Λ)1/2 in A.
*
* @param M input matrix
*
* @param A matrix square root of M
*
* @param lambda the eigenvalues
*
*
*/
public static void PCADecompose (double[][] M, double[][] A,
double[] lambda) {
int d = M.length;
DoubleMatrix2D myM = new DenseDoubleMatrix2D (M);
DoubleMatrix2D myA = new DenseDoubleMatrix2D (d, d);
myA = PCADecompose(myM, lambda);
for(int i = 0; i < d; i++)
for(int j = 0; j < d; j++)
A[i][j] = myA.get(i,j);
}
/**
* Computes the principal components decomposition M =
*
* UΛUt by using the singular value decomposition of matrix
* M. Puts the eigenvalues, which are the diagonal elements of matrix Λ,
* sorted by decreasing size, in vector lambda. Returns matrix
*
* A = U(Λ)1/2.
*
* @param M input matrix
*
* @param lambda the eigenvalues
*
* @return matrix square root of M
*
*/
public static DoubleMatrix2D PCADecompose (DoubleMatrix2D M,
double[] lambda) {
// L'objet SingularValueDecomposition permet de recuperer la matrice
// des valeurs propres en ordre decroissant et celle des vecteurs propres de
// sigma (pour une matrice symétrique et définie-positive seulement).
SingularValueDecomposition sv = new SingularValueDecomposition(M);
// D contient les valeurs propres sur la diagonale
DoubleMatrix2D D = sv.getS ();
for (int i = 0; i < D.rows(); ++i)
lambda[i] = D.getQuick (i, i);
// Calculer la racine carree des valeurs propres
for (int i = 0; i < D.rows(); ++i)
D.setQuick (i, i, Math.sqrt (lambda[i]));
DoubleMatrix2D P = sv.getV();
// Multiplier par la matrice orthogonale (ici P)
return P.zMult (D, null);
}
/**
* Returns matrix M as a string.
* It is displayed in matrix form, with each row on a line.
*
* @return the content of M
*
*/
public static String toString(double[][] M) {
StringBuffer sb = new StringBuffer();
sb.append("{ ");
for(int i = 0; i < M.length; i++) {
sb.append("{ ");
for(int j = 0; j < M[i].length; j++)
sb.append(M[i][j] + " ");
sb.append(" }" + PrintfFormat.NEWLINE + " ");
}
sb.append(" }");
return sb.toString();
}
/**
* Creates a {@link String} containing all the data of
* the DMatrix. The result is displayed in matrix form, with
* each row on a line.
*
* @return the content of the DMatrix
*
*/
public String toString() {
StringBuffer sb = new StringBuffer();
sb.append(" { ");
for(int i = 0; i < r; i++) {
sb.append("{ ");
for(int j = 0; j < c; j++)
sb.append(mat[i][j] + " ");
sb.append(" }" + PrintfFormat.NEWLINE + " ");
}
sb.append("}");
return sb.toString();
}
/**
* Returns the number of rows of the DMatrix.
*
* @return the number of rows
*
*/
public int numRows() {
return r;
}
/**
* Returns the number of columns of the DMatrix.
*
* @return the number of columns
*
*/
public int numColumns() {
return c;
}
/**
* Returns the matrix element in the specified row and column.
*
* @param row the row of the selected element
*
* @param column the column of the selected element
*
* @return the value of the element
* @exception IndexOutOfBoundsException if the selected element would
* be outside the DMatrix
*
*
*/
public double get (int row, int column) {
if (row >= r || column >= c)
throw new IndexOutOfBoundsException();
return mat[row][column];
}
/**
* Sets the value of the element in the specified row and column.
*
* @param row the row of the selected element
*
* @param column the column of the selected element
*
* @param value the new value of the element
*
* @exception IndexOutOfBoundsException if the selected element would
* be outside the DMatrix
*
*
*/
public void set (int row, int column, double value) {
if (row >= r || column >= c)
throw new IndexOutOfBoundsException();
mat[row][column] = value;
}
/**
* Returns the transposed matrix. The rows and columns are
* interchanged.
*
* @return the transposed matrix
*
*/
public DMatrix transpose() {
DMatrix result = new DMatrix(c,r);
for(int i = 0; i < r; i++)
for(int j = 0; j < c; j++)
result.mat[j][i] = mat[i][j];
return result;
}
}
Ads help maintain this website.