org.openimaj.math.matrix.GeneralisedEigenvalueProblem Maven / Gradle / Ivy
Show all versions of cayley Show documentation
/**
* Copyright (c) 2011, The University of Southampton and the individual contributors.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification,
* are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* * Neither the name of the University of Southampton nor the names of its
* contributors may be used to endorse or promote products derived from this
* software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
package org.openimaj.math.matrix;
import gov.nist.math.jama.JamaMatrix;
import no.uib.cipr.matrix.DenseMatrix;
import no.uib.cipr.matrix.DenseVector;
import com.github.fommil.netlib.LAPACK;
import org.netlib.util.intW;
import org.openimaj.math.util.IndependentPair;
import org.openimaj.math.util.ParallelQuicksortDescending;
/**
* Methods for solving the Generalised Eigenvalue Problem: A x = L B x.
*
* Internally the methods in this class use LAPACK to compute the solutions.
*
* @author Jonathon Hare ([email protected])
*/
public final class GeneralisedEigenvalueProblem {
private static int sygvd(int itype, String jobz, String uplo, DenseMatrix A, DenseMatrix B, DenseVector W) {
int info = dsygvd(itype, jobz, uplo, A.numColumns(), A.getData(), A.numRows(), B.getData(), B.numRows(),
W.getData());
if (info == 0) {
return 0;
} else {
if (info < 0) {
throw new IllegalArgumentException("LAPACK ERROR: DSYGVD returned " + info);
}
throw new RuntimeException("LAPACK ERROR: DSYGVD returned " + info);
}
}
private static int dsygvd(int itype, String jobz, String uplo, int n, double[] a, int lda, double[] b, int ldb,
double[] w) {
double[] work = new double[1];
double[] tmp = new double[1];
intW info = new intW(-1);
int lwork;
int[] iwork = new int[1];
int liwork;
// call with info=-1 to determine size of working space
LAPACK.getInstance().dsygvd(itype, jobz, uplo, n, tmp, lda, tmp, ldb, tmp, work, -1, iwork, 0, info);
if (info.val != 0) {
return info.val; // an error occurred
}
// setup working space
lwork = (int) work[0];
work = new double[lwork];
liwork = iwork[0];
iwork = new int[liwork];
// and do the work
LAPACK.getInstance().dsygvd(itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, iwork, liwork, info);
return info.val;
}
/**
* Compute the generalised eigenvalues, L, of the problem A x = L B x. The
* returned eigenvalues are not ordered.
*
* @param A
* symmetric Matrix A; only the upper triangle is used.
* @param B
* symmetric Matrix B; only the upper triangle is used.
* @return the eigenvalues L.
*/
public static DenseVector symmetricGeneralisedEigenvalues(DenseMatrix A, DenseMatrix B) {
if (!A.isSquare() || !B.isSquare()) {
throw new IllegalArgumentException("Input matrices must be square");
}
DenseVector W = new DenseVector(A.numRows());
sygvd(1, "N", "U", A.copy(), B.copy(), W);
return W;
}
/**
* Solve the general problem A x = L B x. The returned eigenvalues are not
* ordered.
*
* @param A
* symmetric matrix A
* @param B
* symmetric matrix B
* @return The eigenvectors x and eigenvalues L.
*/
public static IndependentPair symmetricGeneralisedEigenvectors(DenseMatrix A,
DenseMatrix B) {
if (!A.isSquare() || !B.isSquare()) {
throw new IllegalArgumentException("Input matrices must be square");
}
DenseMatrix vecs = A.copy();
DenseVector W = new DenseVector(A.numRows());
sygvd(1, "V", "U", vecs, B.copy(), W);
return new IndependentPair(vecs, W);
}
/**
* Compute the generalised eigenvalues, L, of the problem A x = L B x. The
* returned eigenvalues are not ordered.
*
* @param A
* symmetric Matrix A; only the upper triangle is used.
* @param B
* symmetric Matrix B; only the upper triangle is used.
* @return the eigenvalues L.
*/
public static double[] symmetricGeneralisedEigenvalues(JamaMatrix A, JamaMatrix B) {
if ((A.getRowDimension() != A.getColumnDimension()) || (B.getRowDimension() != B.getColumnDimension())) {
throw new IllegalArgumentException("Input matrices must be square");
}
DenseVector W = new DenseVector(A.getRowDimension());
sygvd(1, "N", "U", new DenseMatrix(A.getArray()), new DenseMatrix(B.getArray()), W);
return W.getData();
}
/**
* Solve the general problem A x = L B x. The returned eigenvalues are not
* ordered.
*
* @param A
* symmetric matrix A
* @param B
* symmetric matrix B
* @return The eigenvectors x and eigenvalues L.
*/
public static IndependentPair symmetricGeneralisedEigenvectors(JamaMatrix A, JamaMatrix B) {
if ((A.getRowDimension() != A.getColumnDimension()) || (B.getRowDimension() != B.getColumnDimension())) {
throw new IllegalArgumentException("Input matrices must be square");
}
int dim = A.getRowDimension();
DenseMatrix vecs = new DenseMatrix(A.getArray());
DenseVector W = new DenseVector(dim);
sygvd(1, "V", "U", vecs, new DenseMatrix(B.getArray()), W);
JamaMatrix evecs = new JamaMatrix(dim, dim);
double[][] evecsData = evecs.getArray();
double[] vecsData = vecs.getData();
for (int r = 0; r < dim; r++) {
for (int c = 0; c < dim; c++) {
evecsData[r][c] = vecsData[r + c * dim];
}
}
return new IndependentPair(evecs, W.getData());
}
/**
* Solve the general problem A x = L B x. The returned eigenvalues are
* ordered in a decreasing manner.
*
* @param A
* symmetric matrix A
* @param B
* symmetric matrix B
* @return The eigenvectors x and eigenvalues L.
*/
public static IndependentPair symmetricGeneralisedEigenvectorsSorted(JamaMatrix A,
JamaMatrix B) {
if ((A.getRowDimension() != A.getColumnDimension()) || (B.getRowDimension() != B.getColumnDimension())) {
throw new IllegalArgumentException("Input matrices must be square");
}
int dim = A.getRowDimension();
DenseMatrix vecs = new DenseMatrix(A.getArray());
DenseVector W = new DenseVector(dim);
sygvd(1, "V", "U", vecs, new DenseMatrix(B.getArray()), W);
JamaMatrix evecs = new JamaMatrix(dim, dim);
double[][] evecsData = evecs.getArray();
double[] vecsData = vecs.getData();
double[] valsData = W.getData();
int[] indices = range(0, valsData.length - 1);
ParallelQuicksortDescending.sort(valsData, indices);
for (int r = 0; r < dim; r++) {
for (int c = 0; c < dim; c++) {
evecsData[r][c] = vecsData[r + indices[c] * dim];
}
}
return new IndependentPair(evecs, valsData);
}
/**
* Solve the general problem A x = L B x. The returned eigenvalues are
* ordered in a decreasing manner. Only the top numVecs eigenvalues and
* vectors are returned.
*
* @param A
* symmetric matrix A
* @param B
* symmetric matrix B
* @param numVecs
* number of eigenvalues/eigenvectors
* @return The eigenvectors x and eigenvalues L.
*/
public static IndependentPair symmetricGeneralisedEigenvectorsSorted(JamaMatrix A,
JamaMatrix B, int numVecs) {
if ((A.getRowDimension() != A.getColumnDimension()) || (B.getRowDimension() != B.getColumnDimension())) {
throw new IllegalArgumentException("Input matrices must be square");
}
int dim = A.getRowDimension();
DenseMatrix vecs = new DenseMatrix(A.getArray());
DenseVector W = new DenseVector(dim);
sygvd(1, "V", "U", vecs, new DenseMatrix(B.getArray()), W);
JamaMatrix evecs = new JamaMatrix(dim, numVecs);
double[][] evecsData = evecs.getArray();
double[] vecsData = vecs.getData();
double[] valsData = W.getData();
int[] indices = range(0, valsData.length - 1);
ParallelQuicksortDescending.sort(valsData, indices);
for (int r = 0; r < dim; r++) {
for (int c = 0; c < numVecs; c++) {
evecsData[r][c] = vecsData[r + indices[c] * dim];
}
}
return new IndependentPair(evecs, valsData);
}
/**
* Extract a range
*
* @param start
* @param length
* @return [start...length] (inclusive)
*/
private static int[] range(int start, int length) {
int[] range = new int[length - start + 1];
for (int i = start; i <= length; i++) {
range[i - start] = i;
}
return range;
}
}