001package ch.akuhn.matrix.eigenvalues;
002
003import java.util.Random;
004
005import ch.akuhn.matrix.Matrix;
006import ch.akuhn.matrix.SparseMatrix;
007import ch.akuhn.matrix.Vector;
008
009/**
010 * Compute the Singular Value Decomposition.
011 * 
012 * @author Adrian Kuhn
013 */
014public class SingularValues {
015
016        /**
017         * Singular values
018         */
019        public double[] value;
020
021        /**
022         * Left singular vectors
023         */
024        public Vector[] vectorLeft;
025
026        /**
027         * Right singular vectors
028         */
029        public Vector[] vectorRight;
030
031        Matrix A;
032        private int nev;
033
034        /**
035         * Construct with the given matrix and required number of S.V.s
036         * 
037         * @param A
038         * @param nev
039         */
040        public SingularValues(Matrix A, int nev) {
041                this.A = A;
042                this.nev = nev;
043        }
044
045        /**
046         * Perform the decomposition
047         * 
048         * @return this
049         */
050        public SingularValues decompose() {
051                final Eigenvalues eigen = new FewEigenvalues(A.rowCount()) {
052                        @Override
053                        protected Vector callback(Vector v) {
054                                return A.mult(A.transposeMultiply(v));
055                        }
056                }.greatest(nev).run();
057
058                nev = eigen.nev;
059
060                value = new double[nev];
061                vectorLeft = eigen.vector;
062                vectorRight = new Vector[nev];
063                for (int i = 0; i < nev; i++) {
064                        value[i] = Math.sqrt(eigen.value[i]);
065                        vectorRight[i] = A.transposeMultiply(vectorLeft[i]);
066                        vectorRight[i].timesEquals(1.0 / vectorRight[i].norm());
067                }
068                return this;
069        }
070
071        /**
072         * Sample main
073         * 
074         * @param args
075         *            ignored
076         */
077        public static void main(String... args) {
078                final SparseMatrix A = Matrix.sparse(400, 5000);
079                final Random rand = new Random(1);
080                for (int i = 0; i < A.rowCount(); i++) {
081                        for (int j = 0; j < A.columnCount(); j++) {
082                                if (rand.nextDouble() > 0.2)
083                                        continue;
084                                A.put(i, j, rand.nextDouble() * 23);
085                        }
086                }
087
088                final SingularValues singular = new SingularValues(A, 10).decompose();
089                System.out.println(singular.value);
090        }
091}