private static void toR(TridiagonalMatrix tm) { System.out.print("tridiag(c("); for (int i = 0; i < tm.alpha().length; i++) { System.out.print(tm.alpha()[i]); if (i + 1 != tm.alpha().length) { System.out.print(","); } } System.out.print("),c("); for (int i = 0; i < tm.beta().length; i++) { System.out.print(tm.beta()[i]); if (i + 1 != tm.beta().length) { System.out.print(","); } } System.out.println("))"); }
@Override public boolean equals(Object obj) { if (obj == null) { return false; } if (getClass() != obj.getClass()) { return false; } final Solution other = (Solution) obj; if (this.tridiagonal != other.tridiagonal && (this.tridiagonal == null || !this.tridiagonal.equals(other.tridiagonal))) { return false; } if (!Arrays.deepEquals(this.q, other.q)) { return false; } return true; } }
public static <N extends Number> Solution qrSolve(double epsilon, TridiagonalMatrix tridiag, TrivialEigenvalues<N> trivial) { final int n = tridiag.rows(); if (n == 0) { return new Solution(trivial.eigenvalues, new SequenceOfGivens()); } else if (n == 2) { double lambda1, lambda2; lambda1 = lambda2 = tridiag.alpha()[0] + tridiag.alpha()[1]; double t = Math.sqrt((tridiag.alpha()[0] + tridiag.alpha()[1]) * (tridiag.alpha()[0] + tridiag.alpha()[1]) - 4.0 * (tridiag.alpha()[0] * tridiag.alpha()[1] - tridiag.beta()[0] * tridiag.beta()[0])); lambda1 = (lambda1 + t) / 2.0; lambda2 = (lambda2 - t) / 2.0; double u = (lambda1 + tridiag.alpha()[1]) / tridiag.beta()[0]; double c = Math.sqrt(1.0 / (1.0 + u * u)); double s = c / u; if (Math.abs(tridiag.beta()[i]) < epsilon) { if (!qFound) { q++; wilkinsonShift(tridiag.alpha(), tridiag.beta(), p, n - p - q, givensSeq); } else { q++; return new Solution(tridiag.alpha(), givensSeq); } else { double[] eigenvalues = new double[n + trivial.eigenvalues.length]; System.arraycopy(tridiag.alpha(), 0, eigenvalues, 0, n); System.arraycopy(trivial.eigenvalues, 0, eigenvalues, n, trivial.eigenvalues.length); return new Solution(eigenvalues, givensSeq);
@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> Vector<Double> multTransposed(Vector<M> x) { return mult(x); }
public static Solution eigen(VectorFunction<Double, Double> apply, int W, int K, double epsilon) { final LanczosAlgorithm.Solution iLanczos = LanczosAlgorithm.lanczos(apply, randomUnitNormVector(W), K, 1.0); final QRAlgorithm.Solution iQrSolve = QRAlgorithm.qrSolve(epsilon, iLanczos.tridiagonal(), null); final double[][] iEigens = iQrSolve.givensSeq().applyTo(iLanczos.q()); final int[] iOrder = order(iQrSolve.values()); final double[] S = new double[K]; final double[][] V = new double[K][]; for (int i = 0; i < K; i++) { S[i] = Math.sqrt(Math.abs(iLanczos.tridiagonal().alpha()[iOrder[i]])); V[i] = iEigens[iOrder[i]]; } return new Solution(V, null, S); }
beta[i] = Double.parseDouble(ssBeta[i].replaceAll("[\\[\\]\\s]", "")); return new TridiagonalMatrix(alpha, beta);
@Override public int hashCode() { int hash = 5; hash = 71 * hash + (this.tridiagonal != null ? this.tridiagonal.hashCode() : 0); hash = 71 * hash + Arrays.deepHashCode(this.q); return hash; }
@Override public <M extends Number> Vector<Double> mult(Vector<M> x) { return mult(x, Vectors.AS_REALS); }
public static Solution calculateSymmetric(IntIterable corpus, int W, int J, int K, double epsilon) { final LanczosAlgorithm.Solution iLanczos = LanczosAlgorithm.lanczos(new InnerProductMultiplication(corpus, J), randomUnitNormVector(W), K, 1.0); final QRAlgorithm.Solution iQrSolve = QRAlgorithm.qrSolve(epsilon, iLanczos.tridiagonal(), null); final double[][] iEigens = iQrSolve.givensSeq().applyTo(iLanczos.q()); final int[] iOrder = order(iQrSolve.values()); final double[] S = new double[K]; final double[][] V = new double[K][]; for (int i = 0; i < K; i++) { S[i] = Math.sqrt(Math.abs(iLanczos.tridiagonal().alpha()[iOrder[i]])); V[i] = iEigens[iOrder[i]]; } return new Solution(V, null, S); }
assert (K <= n); if (n == 0) { return new Solution(new TridiagonalMatrix(new double[0], new double[0]), new double[0][0], 0.0); } else if (n == 1) { final Vector<Double> unitVector = Vectors.AS_REALS.make(new double[]{1.0}); final double a11 = A.apply(unitVector).doubleValue(0); return new Solution(new TridiagonalMatrix(new double[]{a11}, new double[0]), new double[][]{{1}}, 0.0); } else if (n == 2) { final Vector<Double> unit1Vector = Vectors.AS_REALS.make(new double[]{1.0, 0.0}); final double[] a1 = A.apply(unit1Vector).toDoubleArray(); final double[] a2 = A.apply(unit2Vector).toDoubleArray(); return new Solution(new TridiagonalMatrix(new double[]{a1[0], a2[1]}, new double[]{a1[1]}), new double[][]{ {1, 0}, return new Solution(new TridiagonalMatrix(Arrays.copyOfRange(alpha, 0, K), Arrays.copyOfRange(beta, 1, K)), q, beta[K]);
@Override public Vector<Double> apply(Vector<Double> v) { return mult(v); } };
public static Solution calculate(IntIterable corpus, int W, int J, int K, double epsilon) { final LanczosAlgorithm.Solution oLanczos = LanczosAlgorithm.lanczos(new OuterProductMultiplication(corpus, W), randomUnitNormVector(J), K, 1.0); final QRAlgorithm.Solution oQrSolve = QRAlgorithm.qrSolve(epsilon, oLanczos.tridiagonal(), null); final double[][] oEigens = oQrSolve.givensSeq().applyTo(oLanczos.q()); final int[] oOrder = order(oQrSolve.values()); final LanczosAlgorithm.Solution iLanczos = LanczosAlgorithm.lanczos(new InnerProductMultiplication(corpus, J), randomUnitNormVector(W), K, 1.0); final QRAlgorithm.Solution iQrSolve = QRAlgorithm.qrSolve(epsilon, iLanczos.tridiagonal(), null); final double[][] iEigens = iQrSolve.givensSeq().applyTo(iLanczos.q()); final int[] iOrder = order(iQrSolve.values()); final double[][] V = new double[K][]; final double[][] U = new double[K][]; final double[] S = new double[K]; for (int i = 0; i < K; i++) { S[i] = Math.sqrt(Math.abs((oLanczos.tridiagonal().alpha()[oOrder[i]] + iLanczos.tridiagonal().alpha()[iOrder[i]]) / 2.0)); U[i] = oEigens[oOrder[i]]; V[i] = iEigens[iOrder[i]]; } return new Solution(U, V, S); }
alpha[n - 1] = A2.doubleValue(n - 1, n - 1); return new TridiagonalMatrix(alpha, beta);