/** * Given a symmetric sparse integer matrix compute a tri-diagonalization of * the matrix * * @param A The matrix as a row-wise array of sparse array * @return The diagonal and off-diagonal vector, in such a way that a matrix * B can be constructed such that B[i][i] = r[0][i], B[i][i+1] = r[1][i] and * B[i+1][i] = r[1][i] and B[i][j] = 0 otherwise. This matrix has the same * set of eigenvalues as A. */ public static Solution lanczos(Matrix<Double> A) { assert (A.rows() == A.cols()); assert (A.isSymmetric()); final int n = A.rows(); return lanczos(A.asVectorFunction(), randomUnit(n)); }
@Override public <M extends Number> void sub(Matrix<M> matrix) { for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { data[i][j] -= matrix.doubleValue(i, j); } } }
@Override public Factory<N> factory() { return A.factory(); }
@Override public <M extends Number> void add(Matrix<M> matrix) { assert (matrix.rows() == alpha.length); assert (matrix.cols() == alpha.length); if (matrix instanceof TridiagonalMatrix) { final double[] alpha2 = ((TridiagonalMatrix) matrix).alpha; alpha[i] += matrix.doubleValue(i, i); if (i < alpha.length - 1) { if (matrix.doubleValue(i, i + 1) == matrix.doubleValue(i + 1, i)) { beta[i] += matrix.doubleValue(i, i + 1); } else { throw new RuntimeException("Adding non tridiagonal matrix to tridiagonal matrix!"); if (Math.abs(i - j) > 1 && matrix.doubleValue(i, j) > 1e-6) { throw new RuntimeException("Adding non tridiagonal matrix to tridiagonal matrix!");
assert ((isSymmetric && A.isSymmetric()) || (!isSymmetric && !A.isSymmetric())); final Matrix<N> mat; if (isSymmetric) { mat = A; } else { mat = A.transpose(); final boolean[] trivials = new boolean[mat.rows()]; for (int i = 0; i < mat.rows(); i++) { final Vector<N> row = mat.row(i);
@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); }
/** * Calculate A^TA where A is this matrix * * @param innerProduct The matrix to store the result in * @return */ public void selfInnerProduct(Matrix<N> innerProduct) { assert (innerProduct.cols() == n); assert (innerProduct.rows() == n); for (int k = 0; k < n; k++) { for (Map.Entry<Integer, N> e : arr[k].entrySet()) { final int j = e.getKey(); for (int i = 0; i < n; i++) { final double v = arr[k].doubleValue(i) * e.getValue().doubleValue(); if (v != 0) { innerProduct.add(i, j, v); } } } } }
@Override public <M extends Number> Vector<N> multTransposed(Vector<M> x) { final Vector<N> x2 = A.factory().make(A.rows(), 0.0); for (Map.Entry<Integer, M> e : x.entrySet()) { x2.put(map[e.getKey()], e.getValue().doubleValue()); } return invMap(A.multTransposed(x2), A.factory()); }
return new TrivialEigenvalues<N>(A, new double[0], new boolean[A.rows()]); int[] map = new int[A.rows() - trivColCount]; int trivIdx = 0; int nontrivIdx = 0; for (int i = 0; i < A.rows(); i++) { if (trivialEigenCols[i]) { final Vector<N> row = A.row(i); if (row.size() == 1) { for (Map.Entry<Integer, N> e : row.entrySet()) { if (trivIdx == A.rows()) { final boolean[] allTrue = new boolean[A.rows()]; Arrays.fill(allTrue, false); return new TrivialEigenvalues<N>(null, trivs, allTrue);
public static <N extends Number> TridiagonalMatrix tridiagonalize(Matrix<N> A) { assert (A.isSymmetric()); final int n = A.rows(); final DoubleArrayMatrix A2 = new DoubleArrayMatrix(n, n); A2.add(A);
@Override public boolean isSymmetric() { return A.isSymmetric(); }
@Override public int intValue(int i, int j) { return A.intValue(map[i], map[j]); }
@Override public <M extends Number> void sub(Matrix<M> matrix) { assert (matrix.rows() == alpha.length); assert (matrix.cols() == alpha.length); if (matrix instanceof TridiagonalMatrix) { final double[] alpha2 = ((TridiagonalMatrix) matrix).alpha; alpha[i] -= matrix.doubleValue(i, i); if (i < alpha.length - 1) { if (matrix.doubleValue(i, i + 1) == matrix.doubleValue(i + 1, i)) { beta[i] -= matrix.doubleValue(i, i + 1); } else { throw new RuntimeException("Adding non tridiagonal matrix to tridiagonal matrix!"); if (Math.abs(i - j) > 1 && matrix.doubleValue(i, j) > 1e-6) { throw new RuntimeException("Adding non tridiagonal matrix to tridiagonal matrix!");
/** * Calculate AA^T where A is this matrix * * @param outerProduct The matrix to store the result in * @return */ public void selfOuterProduct(Matrix<N> outerProduct) { assert (outerProduct.cols() == m); assert (outerProduct.rows() == m); for (int i = 0; i < m; i++) { for (Map.Entry<Integer, N> e : arr[i].entrySet()) { for (int j = 0; j < m; j++) { final double v = arr[j].doubleValue(e.getKey()) * e.getValue().doubleValue(); if (v != 0.0) { outerProduct.add(i, j, v); } } } } }
@Override public <M extends Number> Matrix<Double> product(Matrix<M> B) { if (this.cols() != B.rows()) { throw new IllegalArgumentException("Matrix dimensions not suitable for product"); } double[][] res = new double[this.rows()][B.cols()]; for (int i = 0; i < this.rows(); i++) { for (int k = 0; k < B.cols(); k++) { res[i][k] = (i > 0 ? beta[i - 1] * B.doubleValue(i - 1, k) : 0.0) + alpha[i] * B.doubleValue(i, k) + (i + 1 != this.rows() ? beta[i] * B.doubleValue(i + 1, k) : 0.0); } } return new DoubleArrayMatrix(res); } }
@Override public <M extends Number> void add(Matrix<M> matrix) { for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { data[i][j] += matrix.doubleValue(i, j); } } }
@Override public <M extends Number> Matrix<Double> product(Matrix<M> B) { if (this.cols() != B.rows()) { throw new IllegalArgumentException("Matrix dimensions not suitable for product"); } double[][] res = new double[this.rows()][B.cols()]; for (int i = 0; i < this.rows(); i++) { for (int j = 0; j < this.cols(); j++) { for (int k = 0; k < B.cols(); k++) { res[i][k] += data[i][j] * B.doubleValue(j, k); } } } return new DoubleArrayMatrix(res); }
@Override public double doubleValue(int i, int j) { return A.doubleValue(map[i], map[j]); }
@Override public <M extends Number> Matrix<Integer> product(Matrix<M> B) { if (this.cols() != B.rows()) { throw new IllegalArgumentException("Matrix dimensions not suitable for product"); } int[][] res = new int[this.rows()][B.cols()]; for (int i = 0; i < this.rows(); i++) { for (int j = 0; j < this.cols(); j++) { for (int k = 0; k < B.cols(); k++) { res[i][k] += data[i][j] * B.doubleValue(j, k); } } } return new IntArrayMatrix(res); }