public static SparseMatrix<Double> decomp(SparseMatrix<Double> a, boolean complete) { int m = a.rows(); SparseMatrix<Double> l = new SparseMatrix<Double>(m, m, Vectors.AS_SPARSE_REALS); //automatically initialzed to 0's for (int i = 0; i < m; i++) { if (complete) { for (int k = 0; k < (i + 1); k++) { chol(l, k, i, a); } } else { final Vector<Double> l_i = l.row(i); for (int k : l_i.keySet()) { if (k >= (i + 1)) { break; } chol(l, k, i, a); } } } return l; }
@Override public double[][] toDoubleArray() { double[][] d = new double[rows()][cols()]; for (int i = 0; i < rows(); i++) { for (int j = 0; j < rows(); j++) { d[i][j] = doubleValue(i, j); } } return d; } }
@Override public Vector<N> apply(Vector<N> v) { return mult(v); } };
public static void chol(SparseMatrix<Double> l, int k, int i, SparseMatrix<Double> a) throws IllegalArgumentException { double sum = 0; final Vector<Double> l_k = l.row(k); for (int j : l_k.keySet()) { sum += l.doubleValue(i, j) * l_k.doubleValue(j); } double a_ii = a.doubleValue(i, i); if (i == k) { if (a_ii - sum < 0) { throw new IllegalArgumentException("Matrix not positive definite"); } l.set(i, k, Math.sqrt(a_ii - sum)); } else { l.set(i, k, (a.doubleValue(i, k) - sum) / l.doubleValue(k, k)); } }
public static void denseInplaceDecomp(SparseMatrix<Double> a) { int n = a.rows(); int i, j, k; double sum; double[] diagonal = new double[n]; for (i = 0; i < n; i++) { System.err.print("."); for (j = i; j < n; j++) { for (sum = a.doubleValue(i, j), k = i - 1; k >= 0; k--) { sum -= a.doubleValue(i, k) * a.doubleValue(j, k); } if (i == j) { if (sum <= 0.0) { throw new IllegalArgumentException("Matrix not positive definite"); } diagonal[i] = Math.sqrt(sum); } else { a.set(j, i, sum / diagonal[i]); a.set(i, j, 0); } } } for (i = 0; i < n; i++) { a.set(i, i, diagonal[i]); } }
@Override @SuppressWarnings("unchecked") public <M extends Number, O extends Number> Matrix<O> outerProduct(Vector<M> y, Vectors.Factory<O> using) { if (using == Vectors.AS_INTS) { int[][] data2 = new int[n][y.length()]; for (Map.Entry<Integer, Integer> e : entrySet()) { for (int j = 0; j < y.length(); j++) { data2[e.getKey()][j] = e.getValue().intValue() * y.intValue(j); } } return (Matrix<O>) new IntArrayMatrix(data2); } else if (using == Vectors.AS_REALS) { double[][] data2 = new double[n][y.length()]; for (Map.Entry<Integer, Integer> e : entrySet()) { for (int j = 0; j < y.length(); j++) { data2[e.getKey()][j] = y.doubleValue(j) * e.getValue().intValue(); } } return (Matrix<O>) new DoubleArrayMatrix(data2); } else { final SparseMatrix<O> matrix = new SparseMatrix<O>(n, y.length(), using); for (Map.Entry<Integer, Integer> e : entrySet()) { for (Map.Entry<Integer, M> e2 : y.entrySet()) { matrix.set(e.getKey(), e2.getKey(), e2.getValue().doubleValue() * e.getKey().doubleValue()); } } return matrix; } }
@Override public <M extends Number> Matrix<N> product(Matrix<M> B) { if (this.cols() != B.rows()) { throw new IllegalArgumentException("Matrix dimensions not suitable for product"); } if (defaultValue != 0.0 || (B instanceof SparseMatrix && ((SparseMatrix) B).defaultValue != 0.0)) { throw new UnsupportedOperationException(); } Vector<N>[] res = new Vector[this.rows()]; for (int i = 0; i < this.rows(); i++) { res[i] = using.make(B.cols(), 0.0); for (int j : this.arr[i].keySet()) { final Vector<M> r = B.row(j); for (int k : r.keySet()) { res[i].add(k, this.arr[i].doubleValue(j) * B.doubleValue(j, k)); } } } return new SparseMatrix<N>(this.rows(), res, using); }
/** * Find the eigenvalues of sparse matrix by means of a Lanczos-QR method. * * @param A The matrix * @param epsilon The error rate * @return The eigenvalues of the matrix */ public static Solution qrSolve(SparseMatrix<Double> A, double epsilon) { assert (A.rows() == A.cols()); assert (A.isSymmetric()); // It turns out that Lanczos sucks at the easy columns so we remove these // on the first pass final TrivialEigenvalues<Double> trivial = TrivialEigenvalues.find(A, true); if (trivial.nonTrivial == null) { return new Solution(trivial.eigenvalues, new SequenceOfGivens()); } final TridiagonalMatrix tridiag = LanczosAlgorithm.lanczos(trivial.nonTrivial).tridiagonal(); return qrSolve(epsilon, tridiag, trivial); }
public static Vector<Double> solveT(SparseMatrix<Double> a, Vector<Double> b) { assert (a.cols() == b.length()); final int N = b.length(); Vector<Double> y = new SparseRealArray(N); for (int i = N - 1; i >= 0; i--) { double sum = b.doubleValue(i); for (int j = i + 1; j < N; j++) { sum -= a.doubleValue(j, i) * y.doubleValue(j); } y.put(i, sum / a.doubleValue(i, i)); } return y; }
public static <N extends Number> SparseMatrix<N> fromFile(File file, Vectors.Factory<N> using) throws IOException, VectorFormatException { final BufferedReader in = new BufferedReader(new FileReader(file)); String s; final int m = Integer.parseInt(in.readLine()); final int n = Integer.parseInt(in.readLine()); @SuppressWarnings("unchecked") final Vector<N>[] vectors = (Vector<N>[]) new Vector<?>[m]; int idx = 0; while ((s = in.readLine()) != null) { if (idx >= m) { throw new VectorFormatException("Too many lines"); } vectors[idx++] = using.fromString(s, n); } return new SparseMatrix<N>(n, vectors, using); }
@Override public <M extends Number> void add(Matrix<M> matrix) { assert (m == matrix.rows()); assert (n == matrix.cols()); if (matrix instanceof SparseMatrix) { final Vector<M>[] arr2 = ((SparseMatrix<M>) matrix).arr; for (int i = 0; i < m; i++) { if (arr[i] != null && arr2 != null) { arr[i].add(arr2[i]); } else if (arr2 != null) { arr[i] = using.make(n, defaultValue); arr[i].add(arr2[i]); } } } else { for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { row(i).add(j, matrix.doubleValue(i, j)); } } } }
@Override public <M extends Number> Vector<N> multTransposed(Vector<M> x) { assert (m == cols()); double[] result = new double[m]; Arrays.fill(result, m * defaultValue); for (int i = 0; i < arr.length; i++) { for (Map.Entry<Integer, N> e : arr[i].entrySet()) { result[e.getKey().intValue()] += x.doubleValue(i) * e.getValue().doubleValue() - defaultValue; } } return using.make(result); }
@Override @SuppressWarnings("unchecked") public <M extends Number, O extends Number> Matrix<O> outerProduct(Vector<M> y, Vectors.Factory<O> using) { if (using == Vectors.AS_INTS) { int[][] data2 = new int[n][y.length()]; for (Map.Entry<Integer, Double> e : entrySet()) { for (int j = 0; j < y.length(); j++) { data2[e.getKey()][j] = (int) (e.getValue().doubleValue() * y.doubleValue(j)); } } return (Matrix<O>) new IntArrayMatrix(data2); } else if (using == Vectors.AS_REALS) { double[][] data2 = new double[n][y.length()]; for (Map.Entry<Integer, Double> e : entrySet()) { for (int j = 0; j < y.length(); j++) { data2[e.getKey()][j] = y.doubleValue(j) * e.getValue().doubleValue(); } } return (Matrix<O>) new DoubleArrayMatrix(data2); } else { final SparseMatrix<O> matrix = new SparseMatrix<O>(n, y.length(), using); for (Map.Entry<Integer, Double> e : entrySet()) { for (Map.Entry<Integer, M> e2 : y.entrySet()) { matrix.set(e.getKey(), e2.getKey(), e2.getValue().doubleValue() * e.getKey().doubleValue()); } } return matrix; } }
public static Vector<Double> solve(SparseMatrix<Double> a, Vector<Double> b) { assert (a.cols() == b.length()); final int N = b.length(); Vector<Double> y = new SparseRealArray(N); for (int i = 0; i < N; i++) { double sum = b.doubleValue(i); for (int j = 0; j < i; j++) { sum -= a.doubleValue(i, j) * y.doubleValue(j); } y.put(i, sum / a.doubleValue(i, i)); } return y; }
/** * Convert a dense double array to a sparse matrix * * @param arrs The arrays, each element must be non-null and have same * length * @return The sparse matrix */ public static SparseMatrix<Double> fromArray(double[][] arrs) { if (arrs.length == 0) { return new SparseMatrix<Double>(0, 0, AS_SPARSE_REALS); } final SparseMatrix<Double> mat = new SparseMatrix<Double>(arrs.length, arrs[0].length, AS_SPARSE_REALS); for (int i = 0; i < arrs.length; i++) { assert (arrs[i].length == arrs[0].length); for (int j = 0; j < arrs[i].length; j++) { if (arrs[i][j] != 0.0) { if (mat.arr[i] == null) { mat.arr[i] = new SparseRealArray(arrs[i].length); } mat.arr[i].put(j, arrs[i][j]); } } } return mat; }
@Override public <M extends Number> void sub(Matrix<M> matrix) { assert (m == matrix.rows()); assert (n == matrix.cols()); if (matrix instanceof SparseMatrix) { final Vector<M>[] arr2 = ((SparseMatrix<M>) matrix).arr; for (int i = 0; i < m; i++) { if (arr[i] != null && arr2 != null) { arr[i].sub(arr2[i]); } else if (arr2 != null) { arr[i] = using.make(n, defaultValue); arr[i].sub(arr2[i]); } } } else { for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { row(i).sub(j, matrix.doubleValue(i, j)); } } } }
@Override @SuppressWarnings("unchecked") public <M extends Number, O extends Number> Matrix<O> outerProduct(Vector<M> y, Factory<O> using) { if (using == Vectors.AS_INTS) { int[][] data2 = new int[data.length][y.length()]; for (int i = 0; i < data.length; i++) { for (int j = 0; j < y.length(); j++) { data2[i][j] = (int) (data[i] * y.intValue(j)); } } return (Matrix<O>) new IntArrayMatrix(data2); } else if (using == Vectors.AS_REALS) { double[][] data2 = new double[data.length][y.length()]; for (int i = 0; i < data.length; i++) { for (int j = 0; j < y.length(); j++) { data2[i][j] = y.doubleValue(j) * data[i]; } } return (Matrix<O>) new DoubleArrayMatrix(data2); } else { final SparseMatrix<O> matrix = new SparseMatrix<O>(data.length, y.length(), using); for (int i = 0; i < data.length; i++) { for (Map.Entry<Integer, M> e : y.entrySet()) { matrix.set(i, e.getKey(), e.getValue().doubleValue() * data[i]); } } return matrix; } }
/** * Convert a dense double array to a sparse matrix * * @param arrs The arrays, each element must be non-null and have same * length * @return The sparse matrix */ public static SparseMatrix<Integer> fromArray(int[][] arrs) { if (arrs.length == 0) { return new SparseMatrix<Integer>(0, 0, AS_INTS); } final SparseMatrix<Integer> mat = new SparseMatrix<Integer>(arrs.length, arrs[0].length, AS_INTS); for (int i = 0; i < arrs.length; i++) { assert (arrs[i].length == arrs[0].length); for (int j = 0; j < arrs[i].length; j++) { if (arrs[i][j] != 0.0) { if (mat.arr[i] == null) { mat.arr[i] = new SparseIntArray(arrs[i].length); } mat.arr[i].put(j, arrs[i][j]); } } } return mat; }
@Override @SuppressWarnings("unchecked") public <M extends Number, O extends Number> Matrix<O> outerProduct(Vector<M> y, Factory<O> using) { if (using == Vectors.AS_INTS) { int[][] data2 = new int[data.length][y.length()]; for (int i = 0; i < data.length; i++) { for (int j = 0; j < y.length(); j++) { data2[i][j] = data[i] * y.intValue(j); } } return (Matrix<O>) new IntArrayMatrix(data2); } else if (using == Vectors.AS_REALS) { double[][] data2 = new double[data.length][y.length()]; for (int i = 0; i < data.length; i++) { for (int j = 0; j < y.length(); j++) { data2[i][j] = y.doubleValue(j) * data[i]; } } return (Matrix<O>) new DoubleArrayMatrix(data2); } else { final SparseMatrix<O> matrix = new SparseMatrix<O>(data.length, y.length(), using); for (int i = 0; i < data.length; i++) { for (Map.Entry<Integer, M> e : y.entrySet()) { matrix.set(i, e.getKey(), e.getValue().doubleValue() * data[i]); } } return matrix; } }
/** * Generate the transpose of this matrix * * @param using {@link AS_INT_ARRAY} or {@link AS_REAL_ARRAY} as appropriate * @return The transposed matrix */ @Override public SparseMatrix<N> transpose() { final SparseMatrix<N> t = new SparseMatrix<N>(n, m, defaultValue, using); for (int i = 0; i < m; i++) { if (arr[i] != null) { for (Map.Entry<Integer, N> e : arr[i].entrySet()) { final int j = e.getKey(); if (t.arr[j] == null) { t.arr[j] = using.make(m, defaultValue); } t.arr[j].put(i, e.getValue()); } } } return t; }