SingularValues.java

package ch.akuhn.matrix.eigenvalues;

import java.util.Random;

import ch.akuhn.matrix.Matrix;
import ch.akuhn.matrix.SparseMatrix;
import ch.akuhn.matrix.Vector;
import ch.akuhn.util.Out;
import ch.akuhn.util.Stopwatch;

public class SingularValues {

	public double[] value;
	public Vector[] vectorLeft;
	public Vector[] vectorRight;

	Matrix A;
	private int nev;

	public SingularValues(Matrix A, int nev) {
		this.A = A;
		this.nev = nev;
	}

	public SingularValues decompose() {
		Eigenvalues eigen = new FewEigenvalues(A.rowCount()) {
			@Override
			protected Vector callback(Vector v) {
				return A.mult(A.transposeMultiply(v));
			}
		}.greatest(nev).run();
		value = new double[nev];
		vectorLeft = eigen.vector;
		vectorRight = new Vector[nev];
		for (int i = 0; i < nev; i++) {
			value[i] = Math.sqrt(eigen.value[i]);
			vectorRight[i] = A.transposeMultiply(vectorLeft[i]);
			vectorRight[i].times(1.0 / vectorRight[i].norm());
		}
		return this;
	}

	public static void main(String... args) {
		Stopwatch.enter();
		Stopwatch.p();
		SparseMatrix A = Matrix.sparse(4000, 10000);
		Random rand = new Random(1);
		for (int i = 0; i < A.rowCount(); i++) {
			for (int j = 0; j < A.columnCount(); j++) {
				if (rand.nextDouble() > 0.2) continue;
				A.put(i, j, rand.nextDouble() * 23);
			}
		}
		SingularValues singular = new SingularValues(A, 10).decompose();
		Stopwatch.p();
		Out.puts(singular.value);
	}

}