http://git-wip-us.apache.org/repos/asf/ignite/blob/9dccb3db/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/CholeskyDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/CholeskyDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/CholeskyDecomposition.java deleted file mode 100644 index d1a3f51..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/CholeskyDecomposition.java +++ /dev/null @@ -1,309 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.ignite.ml.math.decompositions; - -import org.apache.ignite.ml.math.Destroyable; -import org.apache.ignite.ml.math.Matrix; -import org.apache.ignite.ml.math.Vector; -import org.apache.ignite.ml.math.exceptions.CardinalityException; -import org.apache.ignite.ml.math.exceptions.NonPositiveDefiniteMatrixException; -import org.apache.ignite.ml.math.exceptions.NonSymmetricMatrixException; -import org.apache.ignite.ml.math.util.MatrixUtil; - -import static org.apache.ignite.ml.math.util.MatrixUtil.like; -import static org.apache.ignite.ml.math.util.MatrixUtil.likeVector; - -/** - * Calculates the Cholesky decomposition of a matrix. - * <p> - * This class inspired by class from Apache Common Math with similar name.</p> - * - * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html";>MathWorld</a> - * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition";>Wikipedia</a> - */ -public class CholeskyDecomposition implements Destroyable { - /** - * Default threshold above which off-diagonal elements are considered too different - * and matrix not symmetric. - */ - public static final double DFLT_REL_SYMMETRY_THRESHOLD = 1.0e-15; - - /** - * Default threshold below which diagonal elements are considered null - * and matrix not positive definite. - */ - public static final double DFLT_ABS_POSITIVITY_THRESHOLD = 1.0e-10; - - /** Row-oriented storage for L<sup>T</sup> matrix data. */ - private double[][] lTData; - /** Cached value of L. */ - private Matrix cachedL; - /** Cached value of LT. */ - private Matrix cachedLT; - /** Origin matrix */ - private Matrix origin; - - /** - * Calculates the Cholesky decomposition of the given matrix. - * <p> - * Calling this constructor is equivalent to call {@link #CholeskyDecomposition(Matrix, double, double)} with the - * thresholds set to the default values {@link #DFLT_REL_SYMMETRY_THRESHOLD} and - * {@link #DFLT_ABS_POSITIVITY_THRESHOLD}.</p> - * - * @param mtx the matrix to decompose. - * @throws CardinalityException if matrix is not square. - * @see #CholeskyDecomposition(Matrix, double, double) - * @see #DFLT_REL_SYMMETRY_THRESHOLD - * @see #DFLT_ABS_POSITIVITY_THRESHOLD - */ - public CholeskyDecomposition(final Matrix mtx) { - this(mtx, DFLT_REL_SYMMETRY_THRESHOLD, DFLT_ABS_POSITIVITY_THRESHOLD); - } - - /** - * Calculates the Cholesky decomposition of the given matrix. - * - * @param mtx the matrix to decompose. - * @param relSymmetryThreshold threshold above which off-diagonal elements are considered too different and matrix - * not symmetric - * @param absPositivityThreshold threshold below which diagonal elements are considered null and matrix not positive - * definite - * @see #CholeskyDecomposition(Matrix) - * @see #DFLT_REL_SYMMETRY_THRESHOLD - * @see #DFLT_ABS_POSITIVITY_THRESHOLD - */ - public CholeskyDecomposition(final Matrix mtx, final double relSymmetryThreshold, - final double absPositivityThreshold) { - assert mtx != null; - - if (mtx.columnSize() != mtx.rowSize()) - throw new CardinalityException(mtx.rowSize(), mtx.columnSize()); - - origin = mtx; - - final int order = mtx.rowSize(); - - lTData = toDoubleArr(mtx); - cachedL = null; - cachedLT = null; - - // Check the matrix before transformation. - for (int i = 0; i < order; ++i) { - final double[] lI = lTData[i]; - - // Check off-diagonal elements (and reset them to 0). - for (int j = i + 1; j < order; ++j) { - final double[] lJ = lTData[j]; - - final double lIJ = lI[j]; - final double lJI = lJ[i]; - - final double maxDelta = relSymmetryThreshold * Math.max(Math.abs(lIJ), Math.abs(lJI)); - - if (Math.abs(lIJ - lJI) > maxDelta) - throw new NonSymmetricMatrixException(i, j, relSymmetryThreshold); - - lJ[i] = 0; - } - } - - // Transform the matrix. - for (int i = 0; i < order; ++i) { - final double[] ltI = lTData[i]; - - // Check diagonal element. - if (ltI[i] <= absPositivityThreshold) - throw new NonPositiveDefiniteMatrixException(ltI[i], i, absPositivityThreshold); - - ltI[i] = Math.sqrt(ltI[i]); - final double inverse = 1.0 / ltI[i]; - - for (int q = order - 1; q > i; --q) { - ltI[q] *= inverse; - final double[] ltQ = lTData[q]; - - for (int p = q; p < order; ++p) - ltQ[p] -= ltI[q] * ltI[p]; - } - } - } - - /** {@inheritDoc} */ - @Override public void destroy() { - if (cachedL != null) - cachedL.destroy(); - if (cachedLT != null) - cachedLT.destroy(); - } - - /** - * Returns the matrix L of the decomposition. - * <p>L is an lower-triangular matrix</p> - * - * @return the L matrix - */ - public Matrix getL() { - if (cachedL == null) - cachedL = getLT().transpose(); - - return cachedL; - } - - /** - * Returns the transpose of the matrix L of the decomposition. - * <p>L<sup>T</sup> is an upper-triangular matrix</p> - * - * @return the transpose of the matrix L of the decomposition - */ - public Matrix getLT() { - - if (cachedLT == null) { - Matrix like = like(origin, origin.rowSize(), origin.columnSize()); - like.assign(lTData); - - cachedLT = like; - } - - // return the cached matrix - return cachedLT; - } - - /** - * Return the determinant of the matrix - * - * @return determinant of the matrix - */ - public double getDeterminant() { - double determinant = 1.0; - - for (int i = 0; i < lTData.length; ++i) { - double lTii = lTData[i][i]; - determinant *= lTii * lTii; - } - - return determinant; - } - - /** - * Solve the linear equation A × X = B for matrices A. - * - * @param b right-hand side of the equation A × X = B - * @return a vector X that minimizes the two norm of A × X - B - * @throws CardinalityException if the vectors dimensions do not match - */ - public Vector solve(final Vector b) { - final int m = lTData.length; - - if (b.size() != m) - throw new CardinalityException(b.size(), m); - - final double[] x = b.getStorage().data(); - - // Solve LY = b - for (int j = 0; j < m; j++) { - final double[] lJ = lTData[j]; - - x[j] /= lJ[j]; - - final double xJ = x[j]; - - for (int i = j + 1; i < m; i++) - x[i] -= xJ * lJ[i]; - } - - // Solve LTX = Y - for (int j = m - 1; j >= 0; j--) { - x[j] /= lTData[j][j]; - - final double xJ = x[j]; - - for (int i = 0; i < j; i++) - x[i] -= xJ * lTData[i][j]; - } - - return likeVector(origin, m).assign(x); - } - - /** - * Solve the linear equation A × X = B for matrices A. - * - * @param b right-hand side of the equation A × X = B - * @return a matrix X that minimizes the two norm of A × X - B - * @throws CardinalityException if the matrices dimensions do not match - */ - public Matrix solve(final Matrix b) { - final int m = lTData.length; - - if (b.rowSize() != m) - throw new CardinalityException(b.rowSize(), m); - - final int nColB = b.columnSize(); - final double[][] x = MatrixUtil.unflatten(b.getStorage().data(), b.columnSize(), b.getStorage().storageMode()); - - // Solve LY = b - for (int j = 0; j < m; j++) { - final double[] lJ = lTData[j]; - final double lJJ = lJ[j]; - final double[] xJ = x[j]; - - for (int k = 0; k < nColB; ++k) - xJ[k] /= lJJ; - - for (int i = j + 1; i < m; i++) { - final double[] xI = x[i]; - final double lJI = lJ[i]; - - for (int k = 0; k < nColB; ++k) - xI[k] -= xJ[k] * lJI; - } - } - - // Solve LTX = Y - for (int j = m - 1; j >= 0; j--) { - final double lJJ = lTData[j][j]; - final double[] xJ = x[j]; - - for (int k = 0; k < nColB; ++k) - xJ[k] /= lJJ; - - for (int i = 0; i < j; i++) { - final double[] xI = x[i]; - final double lIJ = lTData[i][j]; - - for (int k = 0; k < nColB; ++k) - xI[k] -= xJ[k] * lIJ; - } - } - - return like(origin, m, b.columnSize()).assign(x); - } - - /** */ - private double[][] toDoubleArr(Matrix mtx) { - if (mtx.isArrayBased()) - return MatrixUtil.unflatten(mtx.getStorage().data(), mtx.columnSize(), mtx.getStorage().storageMode()); - - double[][] res = new double[mtx.rowSize()][mtx.columnSize()]; - - for (int row = 0; row < mtx.rowSize(); row++) - for (int col = 0; col < mtx.columnSize(); col++) - res[row][col] = mtx.get(row, col); - - return res; - } -}
http://git-wip-us.apache.org/repos/asf/ignite/blob/9dccb3db/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/DecompositionSupport.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/DecompositionSupport.java b/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/DecompositionSupport.java deleted file mode 100644 index 20d6e79..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/DecompositionSupport.java +++ /dev/null @@ -1,105 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.ignite.ml.math.decompositions; - -import org.apache.ignite.ml.math.Destroyable; -import org.apache.ignite.ml.math.Matrix; -import org.apache.ignite.ml.math.Vector; -import org.apache.ignite.ml.math.impls.matrix.CacheMatrix; -import org.apache.ignite.ml.math.impls.matrix.DenseLocalOnHeapMatrix; -import org.apache.ignite.ml.math.impls.matrix.PivotedMatrixView; -import org.apache.ignite.ml.math.impls.matrix.RandomMatrix; -import org.apache.ignite.ml.math.impls.vector.DenseLocalOnHeapVector; - -/** - * Helper methods to support decomposition of matrix types having some functionality limited. - */ -public abstract class DecompositionSupport implements Destroyable { - /** - * Create the like matrix with read-only matrices support. - * - * @param matrix Matrix for like. - * @return Like matrix. - */ - protected Matrix like(Matrix matrix) { - if (isCopyLikeSupport(matrix)) - return new DenseLocalOnHeapMatrix(matrix.rowSize(), matrix.columnSize()); - else - return matrix.like(matrix.rowSize(), matrix.columnSize()); - } - - /** - * Create the like matrix with specified size with read-only matrices support. - * - * @param matrix Matrix for like. - * @return Like matrix. - */ - protected Matrix like(Matrix matrix, int rows, int cols) { - if (isCopyLikeSupport(matrix)) - return new DenseLocalOnHeapMatrix(rows, cols); - else - return matrix.like(rows, cols); - } - - /** - * Create the like vector with read-only matrices support. - * - * @param matrix Matrix for like. - * @param crd Cardinality of the vector. - * @return Like vector. - */ - protected Vector likeVector(Matrix matrix, int crd) { - if (isCopyLikeSupport(matrix)) - return new DenseLocalOnHeapVector(crd); - else - return matrix.likeVector(crd); - } - - /** - * Create the like vector with read-only matrices support. - * - * @param matrix Matrix for like. - * @return Like vector. - */ - protected Vector likeVector(Matrix matrix) { - return likeVector(matrix, matrix.rowSize()); - } - - /** - * Create the copy of matrix with read-only matrices support. - * - * @param matrix Matrix for copy. - * @return Copy. - */ - protected Matrix copy(Matrix matrix) { - if (isCopyLikeSupport(matrix)) { - DenseLocalOnHeapMatrix cp = new DenseLocalOnHeapMatrix(matrix.rowSize(), matrix.columnSize()); - - cp.assign(matrix); - - return cp; - } - else - return matrix.copy(); - } - - /** */ - private boolean isCopyLikeSupport(Matrix matrix) { - return matrix instanceof RandomMatrix || matrix instanceof PivotedMatrixView || matrix instanceof CacheMatrix; - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/9dccb3db/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/EigenDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/EigenDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/EigenDecomposition.java deleted file mode 100644 index a5c92e6..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/EigenDecomposition.java +++ /dev/null @@ -1,936 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.ignite.ml.math.decompositions; - -import org.apache.ignite.ml.math.Destroyable; -import org.apache.ignite.ml.math.Matrix; -import org.apache.ignite.ml.math.Vector; -import org.apache.ignite.ml.math.functions.Functions; - -import static org.apache.ignite.ml.math.util.MatrixUtil.like; -import static org.apache.ignite.ml.math.util.MatrixUtil.likeVector; - -/** - * This class provides EigenDecomposition of given matrix. The class is based on - * class with similar name from <a href="http://mahout.apache.org/";>Apache Mahout</a> library. - * - * @see <a href=http://mathworld.wolfram.com/EigenDecomposition.html>MathWorld</a> - */ -public class EigenDecomposition implements Destroyable { - /** Row and column dimension (square matrix). */ - private final int n; - - /** Array for internal storage of eigen vectors. */ - private final Matrix v; - - /** Array for internal storage of eigenvalues. */ - private final Vector d; - /** Array for internal storage of eigenvalues. */ - private final Vector e; - - /** - * @param matrix Matrix to decompose. - */ - public EigenDecomposition(Matrix matrix) { - this(matrix, isSymmetric(matrix)); - } - - /** - * @param matrix Matrix to decompose. - * @param isSymmetric {@code true} if matrix passes symmetry check, {@code false otherwise}. - */ - public EigenDecomposition(Matrix matrix, boolean isSymmetric) { - n = matrix.columnSize(); - - d = likeVector(matrix); - e = likeVector(matrix); - v = like(matrix); - - if (isSymmetric) { - v.assign(matrix); - - // Tridiagonalize. - tred2(); - - // Diagonalize. - tql2(); - - } - else - // Reduce to Hessenberg form. - // Reduce Hessenberg to real Schur form. - hqr2(orthes(matrix)); - } - - /** - * Return the eigen vector matrix. - * - * @return V - */ - public Matrix getV() { - return like(v).assign(v); - } - - /** - * Return the real parts of the eigenvalues - */ - public Vector getRealEigenValues() { - return d; - } - - /** - * Return the imaginary parts of the eigenvalues. - * - * @return Vector of imaginary parts. - */ - public Vector getImagEigenvalues() { - return e; - } - - /** - * Return the block diagonal eigenvalue matrix - * - * @return D - */ - public Matrix getD() { - Matrix res = like(v, d.size(), d.size()); - res.assign(0); - res.viewDiagonal().assign(d); - for (int i = 0; i < n; i++) { - double v = e.getX(i); - if (v > 0) - res.setX(i, i + 1, v); - else if (v < 0) - res.setX(i, i - 1, v); - } - return res; - } - - /** - * Destroys decomposition components and other internal components of decomposition. - */ - @Override public void destroy() { - e.destroy(); - v.destroy(); - d.destroy(); - } - - /** */ - private void tred2() { - // This is derived from the Algol procedures tred2 by - // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for - // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - d.assign(v.viewColumn(n - 1)); - - // Householder reduction to tridiagonal form. - - for (int i = n - 1; i > 0; i--) { - - // Scale to avoid under/overflow. - double scale = d.viewPart(0, i).kNorm(1); - double h = 0.0; - - if (scale == 0.0) { - e.setX(i, d.getX(i - 1)); - for (int j = 0; j < i; j++) { - d.setX(j, v.getX(i - 1, j)); - v.setX(i, j, 0.0); - v.setX(j, i, 0.0); - } - } - else { - - // Generate Householder vector. - - for (int k = 0; k < i; k++) { - d.setX(k, d.getX(k) / scale); - h += d.getX(k) * d.getX(k); - } - - double f = d.getX(i - 1); - double g = Math.sqrt(h); - - if (f > 0) - g = -g; - - e.setX(i, scale * g); - h -= f * g; - d.setX(i - 1, f - g); - - for (int j = 0; j < i; j++) - e.setX(j, 0.0); - - // Apply similarity transformation to remaining columns. - - for (int j = 0; j < i; j++) { - f = d.getX(j); - v.setX(j, i, f); - g = e.getX(j) + v.getX(j, j) * f; - - for (int k = j + 1; k <= i - 1; k++) { - g += v.getX(k, j) * d.getX(k); - e.setX(k, e.getX(k) + v.getX(k, j) * f); - } - - e.setX(j, g); - } - - f = 0.0; - - for (int j = 0; j < i; j++) { - e.setX(j, e.getX(j) / h); - f += e.getX(j) * d.getX(j); - } - - double hh = f / (h + h); - - for (int j = 0; j < i; j++) - e.setX(j, e.getX(j) - hh * d.getX(j)); - - for (int j = 0; j < i; j++) { - f = d.getX(j); - g = e.getX(j); - - for (int k = j; k <= i - 1; k++) - v.setX(k, j, v.getX(k, j) - (f * e.getX(k) + g * d.getX(k))); - - d.setX(j, v.getX(i - 1, j)); - v.setX(i, j, 0.0); - } - } - - d.setX(i, h); - } - } - - /** */ - private Matrix orthes(Matrix matrix) { - // Working storage for nonsymmetric algorithm. - Vector ort = likeVector(matrix); - Matrix hessenBerg = like(matrix).assign(matrix); - - // This is derived from the Algol procedures orthes and ortran, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutines in EISPACK. - - int low = 0; - int high = n - 1; - - for (int m = low + 1; m <= high - 1; m++) { - - // Scale column. - - Vector hCol = hessenBerg.viewColumn(m - 1).viewPart(m, high - m + 1); - double scale = hCol.kNorm(1); - - if (scale != 0.0) { - // Compute Householder transformation. - ort.viewPart(m, high - m + 1).map(hCol, Functions.plusMult(1 / scale)); - double h = ort.viewPart(m, high - m + 1).getLengthSquared(); - - double g = Math.sqrt(h); - - if (ort.getX(m) > 0) - g = -g; - - h -= ort.getX(m) * g; - ort.setX(m, ort.getX(m) - g); - - // Apply Householder similarity transformation - // H = (I-u*u'/h)*H*(I-u*u')/h) - - Vector ortPiece = ort.viewPart(m, high - m + 1); - - for (int j = m; j < n; j++) { - double f = ortPiece.dot(hessenBerg.viewColumn(j).viewPart(m, high - m + 1)) / h; - hessenBerg.viewColumn(j).viewPart(m, high - m + 1).map(ortPiece, Functions.plusMult(-f)); - } - - for (int i = 0; i <= high; i++) { - double f = ortPiece.dot(hessenBerg.viewRow(i).viewPart(m, high - m + 1)) / h; - hessenBerg.viewRow(i).viewPart(m, high - m + 1).map(ortPiece, Functions.plusMult(-f)); - } - - ort.setX(m, scale * ort.getX(m)); - hessenBerg.setX(m, m - 1, scale * g); - } - } - - // Accumulate transformations (Algol's ortran). - - v.assign(0); - v.viewDiagonal().assign(1); - - for (int m = high - 1; m >= low + 1; m--) { - if (hessenBerg.getX(m, m - 1) != 0.0) { - ort.viewPart(m + 1, high - m).assign(hessenBerg.viewColumn(m - 1).viewPart(m + 1, high - m)); - - for (int j = m; j <= high; j++) { - double g = ort.viewPart(m, high - m + 1).dot(v.viewColumn(j).viewPart(m, high - m + 1)); - - // Double division avoids possible underflow - g = g / ort.getX(m) / hessenBerg.getX(m, m - 1); - v.viewColumn(j).viewPart(m, high - m + 1).map(ort.viewPart(m, high - m + 1), Functions.plusMult(g)); - } - } - } - - return hessenBerg; - } - - /** - * Symmetric tridiagonal QL algorithm. - */ - private void tql2() { - // This is derived from the Algol procedures tql2, by - // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for - // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - e.viewPart(0, n - 1).assign(e.viewPart(1, n - 1)); - e.setX(n - 1, 0.0); - - double f = 0.0; - double tst1 = 0.0; - double eps = Math.pow(2.0, -52.0); - - for (int l = 0; l < n; l++) { - // Find small subdiagonal element. - - tst1 = Math.max(tst1, Math.abs(d.getX(l)) + Math.abs(e.getX(l))); - int m = l; - - while (m < n) { - if (Math.abs(e.getX(m)) <= eps * tst1) - break; - - m++; - } - - // If m == l, d.getX(l) is an eigenvalue, - // otherwise, iterate. - - if (m > l) { - do { - // Compute implicit shift - - double g = d.getX(l); - double p = (d.getX(l + 1) - g) / (2.0 * e.getX(l)); - double r = Math.hypot(p, 1.0); - - if (p < 0) - r = -r; - - d.setX(l, e.getX(l) / (p + r)); - d.setX(l + 1, e.getX(l) * (p + r)); - double dl1 = d.getX(l + 1); - double h = g - d.getX(l); - - for (int i = l + 2; i < n; i++) - d.setX(i, d.getX(i) - h); - - f += h; - - // Implicit QL transformation. - - p = d.getX(m); - double c = 1.0; - double c2 = c; - double c3 = c; - double el1 = e.getX(l + 1); - double s = 0.0; - double s2 = 0.0; - - for (int i = m - 1; i >= l; i--) { - c3 = c2; - c2 = c; - s2 = s; - g = c * e.getX(i); - h = c * p; - r = Math.hypot(p, e.getX(i)); - e.setX(i + 1, s * r); - s = e.getX(i) / r; - c = p / r; - p = c * d.getX(i) - s * g; - d.setX(i + 1, h + s * (c * g + s * d.getX(i))); - - // Accumulate transformation. - - for (int k = 0; k < n; k++) { - h = v.getX(k, i + 1); - v.setX(k, i + 1, s * v.getX(k, i) + c * h); - v.setX(k, i, c * v.getX(k, i) - s * h); - } - } - - p = -s * s2 * c3 * el1 * e.getX(l) / dl1; - e.setX(l, s * p); - d.setX(l, c * p); - - // Check for convergence. - - } - while (Math.abs(e.getX(l)) > eps * tst1); - } - - d.setX(l, d.getX(l) + f); - e.setX(l, 0.0); - } - - // Sort eigenvalues and corresponding vectors. - - for (int i = 0; i < n - 1; i++) { - int k = i; - double p = d.getX(i); - - for (int j = i + 1; j < n; j++) - if (d.getX(j) > p) { - k = j; - p = d.getX(j); - } - - if (k != i) { - d.setX(k, d.getX(i)); - d.setX(i, p); - - for (int j = 0; j < n; j++) { - p = v.getX(j, i); - v.setX(j, i, v.getX(j, k)); - v.setX(j, k, p); - } - } - } - } - - /** */ - private void hqr2(Matrix h) { - // This is derived from the Algol procedure hqr2, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. - - // Initialize - - int nn = this.n; - int n = nn - 1; - int low = 0; - int high = nn - 1; - double eps = Math.pow(2.0, -52.0); - double exshift = 0.0; - double p = 0; - double q = 0; - double r = 0; - double s = 0; - double z = 0; - double w; - double x; - double y; - - // Store roots isolated by balanc and compute matrix norm - - double norm = h.foldMap(Functions.PLUS, Functions.ABS, 0.0d); - - // Outer loop over eigenvalue index - - int iter = 0; - while (n >= low) { - // Look for single small sub-diagonal element - int l = n; - - while (l > low) { - s = Math.abs(h.getX(l - 1, l - 1)) + Math.abs(h.getX(l, l)); - - if (s == 0.0) - s = norm; - - if (Math.abs(h.getX(l, l - 1)) < eps * s) - break; - - l--; - } - - // Check for convergence - - if (l == n) { - // One root found - h.setX(n, n, h.getX(n, n) + exshift); - d.setX(n, h.getX(n, n)); - e.setX(n, 0.0); - n--; - iter = 0; - } - else if (l == n - 1) { - // Two roots found - w = h.getX(n, n - 1) * h.getX(n - 1, n); - p = (h.getX(n - 1, n - 1) - h.getX(n, n)) / 2.0; - q = p * p + w; - z = Math.sqrt(Math.abs(q)); - h.setX(n, n, h.getX(n, n) + exshift); - h.setX(n - 1, n - 1, h.getX(n - 1, n - 1) + exshift); - x = h.getX(n, n); - - // Real pair - if (q >= 0) { - if (p >= 0) - z = p + z; - else - z = p - z; - - d.setX(n - 1, x + z); - d.setX(n, d.getX(n - 1)); - - if (z != 0.0) - d.setX(n, x - w / z); - - e.setX(n - 1, 0.0); - e.setX(n, 0.0); - x = h.getX(n, n - 1); - s = Math.abs(x) + Math.abs(z); - p = x / s; - q = z / s; - r = Math.sqrt(p * p + q * q); - p /= r; - q /= r; - - // Row modification - - for (int j = n - 1; j < nn; j++) { - z = h.getX(n - 1, j); - h.setX(n - 1, j, q * z + p * h.getX(n, j)); - h.setX(n, j, q * h.getX(n, j) - p * z); - } - - // Column modification - - for (int i = 0; i <= n; i++) { - z = h.getX(i, n - 1); - h.setX(i, n - 1, q * z + p * h.getX(i, n)); - h.setX(i, n, q * h.getX(i, n) - p * z); - } - - // Accumulate transformations - - for (int i = low; i <= high; i++) { - z = v.getX(i, n - 1); - v.setX(i, n - 1, q * z + p * v.getX(i, n)); - v.setX(i, n, q * v.getX(i, n) - p * z); - } - - // Complex pair - - } - else { - d.setX(n - 1, x + p); - d.setX(n, x + p); - e.setX(n - 1, z); - e.setX(n, -z); - } - - n -= 2; - iter = 0; - - // No convergence yet - - } - else { - // Form shift - x = h.getX(n, n); - y = 0.0; - w = 0.0; - - if (l < n) { - y = h.getX(n - 1, n - 1); - w = h.getX(n, n - 1) * h.getX(n - 1, n); - } - - // Wilkinson's original ad hoc shift - - if (iter == 10) { - exshift += x; - - for (int i = low; i <= n; i++) - h.setX(i, i, x); - - s = Math.abs(h.getX(n, n - 1)) + Math.abs(h.getX(n - 1, n - 2)); - x = y = 0.75 * s; - w = -0.4375 * s * s; - } - - // MATLAB's new ad hoc shift - - if (iter == 30) { - s = (y - x) / 2.0; - s = s * s + w; - - if (s > 0) { - s = Math.sqrt(s); - - if (y < x) - s = -s; - - s = x - w / ((y - x) / 2.0 + s); - - for (int i = low; i <= n; i++) - h.setX(i, i, h.getX(i, i) - s); - - exshift += s; - x = y = w = 0.964; - } - } - - iter++; // (Could check iteration count here.) - - // Look for two consecutive small sub-diagonal elements - - int m = n - 2; - - while (m >= l) { - z = h.getX(m, m); - r = x - z; - s = y - z; - p = (r * s - w) / h.getX(m + 1, m) + h.getX(m, m + 1); - q = h.getX(m + 1, m + 1) - z - r - s; - r = h.getX(m + 2, m + 1); - s = Math.abs(p) + Math.abs(q) + Math.abs(r); - p /= s; - q /= s; - r /= s; - - if (m == l) - break; - - double hmag = Math.abs(h.getX(m - 1, m - 1)) + Math.abs(h.getX(m + 1, m + 1)); - double threshold = eps * Math.abs(p) * (Math.abs(z) + hmag); - - if (Math.abs(h.getX(m, m - 1)) * (Math.abs(q) + Math.abs(r)) < threshold) - break; - - m--; - } - - for (int i = m + 2; i <= n; i++) { - h.setX(i, i - 2, 0.0); - - if (i > m + 2) - h.setX(i, i - 3, 0.0); - } - - // Double QR step involving rows l:n and columns m:n - - for (int k = m; k <= n - 1; k++) { - boolean notlast = k != n - 1; - - if (k != m) { - p = h.getX(k, k - 1); - q = h.getX(k + 1, k - 1); - r = notlast ? h.getX(k + 2, k - 1) : 0.0; - x = Math.abs(p) + Math.abs(q) + Math.abs(r); - if (x != 0.0) { - p /= x; - q /= x; - r /= x; - } - } - - if (x == 0.0) - break; - - s = Math.sqrt(p * p + q * q + r * r); - - if (p < 0) - s = -s; - - if (s != 0) { - if (k != m) - h.setX(k, k - 1, -s * x); - else if (l != m) - h.setX(k, k - 1, -h.getX(k, k - 1)); - - p += s; - x = p / s; - y = q / s; - z = r / s; - q /= p; - r /= p; - - // Row modification - - for (int j = k; j < nn; j++) { - p = h.getX(k, j) + q * h.getX(k + 1, j); - - if (notlast) { - p += r * h.getX(k + 2, j); - h.setX(k + 2, j, h.getX(k + 2, j) - p * z); - } - - h.setX(k, j, h.getX(k, j) - p * x); - h.setX(k + 1, j, h.getX(k + 1, j) - p * y); - } - - // Column modification - - for (int i = 0; i <= Math.min(n, k + 3); i++) { - p = x * h.getX(i, k) + y * h.getX(i, k + 1); - - if (notlast) { - p += z * h.getX(i, k + 2); - h.setX(i, k + 2, h.getX(i, k + 2) - p * r); - } - - h.setX(i, k, h.getX(i, k) - p); - h.setX(i, k + 1, h.getX(i, k + 1) - p * q); - } - - // Accumulate transformations - - for (int i = low; i <= high; i++) { - p = x * v.getX(i, k) + y * v.getX(i, k + 1); - - if (notlast) { - p += z * v.getX(i, k + 2); - v.setX(i, k + 2, v.getX(i, k + 2) - p * r); - } - - v.setX(i, k, v.getX(i, k) - p); - v.setX(i, k + 1, v.getX(i, k + 1) - p * q); - } - } // (s != 0) - } // k loop - } // check convergence - } // while (n >= low) - - // Back substitute to find vectors of upper triangular form - - if (norm == 0.0) - return; - - for (n = nn - 1; n >= 0; n--) { - p = d.getX(n); - q = e.getX(n); - - // Real vector - - double t; - - if (q == 0) { - int l = n; - h.setX(n, n, 1.0); - - for (int i = n - 1; i >= 0; i--) { - w = h.getX(i, i) - p; - r = 0.0; - - for (int j = l; j <= n; j++) - r += h.getX(i, j) * h.getX(j, n); - - if (e.getX(i) < 0.0) { - z = w; - s = r; - } - else { - l = i; - - if (e.getX(i) == 0.0) { - if (w == 0.0) - h.setX(i, n, -r / (eps * norm)); - else - h.setX(i, n, -r / w); - - // Solve real equations - - } - else { - x = h.getX(i, i + 1); - y = h.getX(i + 1, i); - q = (d.getX(i) - p) * (d.getX(i) - p) + e.getX(i) * e.getX(i); - t = (x * s - z * r) / q; - h.setX(i, n, t); - - if (Math.abs(x) > Math.abs(z)) - h.setX(i + 1, n, (-r - w * t) / x); - else - h.setX(i + 1, n, (-s - y * t) / z); - } - - // Overflow control - - t = Math.abs(h.getX(i, n)); - - if (eps * t * t > 1) { - for (int j = i; j <= n; j++) - h.setX(j, n, h.getX(j, n) / t); - } - } - } - - // Complex vector - - } - else if (q < 0) { - int l = n - 1; - - // Last vector component imaginary so matrix is triangular - - if (Math.abs(h.getX(n, n - 1)) > Math.abs(h.getX(n - 1, n))) { - h.setX(n - 1, n - 1, q / h.getX(n, n - 1)); - h.setX(n - 1, n, -(h.getX(n, n) - p) / h.getX(n, n - 1)); - } - else { - cdiv(0.0, -h.getX(n - 1, n), h.getX(n - 1, n - 1) - p, q); - h.setX(n - 1, n - 1, cdivr); - h.setX(n - 1, n, cdivi); - } - - h.setX(n, n - 1, 0.0); - h.setX(n, n, 1.0); - - for (int i = n - 2; i >= 0; i--) { - double ra = 0.0; - double sa = 0.0; - - for (int j = l; j <= n; j++) { - ra += h.getX(i, j) * h.getX(j, n - 1); - sa += h.getX(i, j) * h.getX(j, n); - } - - w = h.getX(i, i) - p; - - if (e.getX(i) < 0.0) { - z = w; - r = ra; - s = sa; - } - else { - l = i; - - if (e.getX(i) == 0) { - cdiv(-ra, -sa, w, q); - h.setX(i, n - 1, cdivr); - h.setX(i, n, cdivi); - } - else { - - // Solve complex equations - - x = h.getX(i, i + 1); - y = h.getX(i + 1, i); - - double vr = (d.getX(i) - p) * (d.getX(i) - p) + e.getX(i) * e.getX(i) - q * q; - double vi = (d.getX(i) - p) * 2.0 * q; - - if (vr == 0.0 && vi == 0.0) { - double hmag = Math.abs(x) + Math.abs(y); - vr = eps * norm * (Math.abs(w) + Math.abs(q) + hmag + Math.abs(z)); - } - - cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi); - - h.setX(i, n - 1, cdivr); - h.setX(i, n, cdivi); - - if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) { - h.setX(i + 1, n - 1, (-ra - w * h.getX(i, n - 1) + q * h.getX(i, n)) / x); - h.setX(i + 1, n, (-sa - w * h.getX(i, n) - q * h.getX(i, n - 1)) / x); - } - else { - cdiv(-r - y * h.getX(i, n - 1), -s - y * h.getX(i, n), z, q); - - h.setX(i + 1, n - 1, cdivr); - h.setX(i + 1, n, cdivi); - } - } - - // Overflow control - - t = Math.max(Math.abs(h.getX(i, n - 1)), Math.abs(h.getX(i, n))); - - if (eps * t * t > 1) - for (int j = i; j <= n; j++) { - h.setX(j, n - 1, h.getX(j, n - 1) / t); - h.setX(j, n, h.getX(j, n) / t); - } - } - } - } - } - - // Vectors of isolated roots - - for (int i = 0; i < nn; i++) - if (i < low || i > high) { - for (int j = i; j < nn; j++) - v.setX(i, j, h.getX(i, j)); - } - - // Back transformation to get eigen vectors of original matrix - - for (int j = nn - 1; j >= low; j--) - for (int i = low; i <= high; i++) { - z = 0.0; - - for (int k = low; k <= Math.min(j, high); k++) - z += v.getX(i, k) * h.getX(k, j); - - v.setX(i, j, z); - } - } - - /** */ - private static boolean isSymmetric(Matrix matrix) { - int cols = matrix.columnSize(); - int rows = matrix.rowSize(); - - if (cols != rows) - return false; - - for (int i = 0; i < cols; i++) - for (int j = 0; j < rows; j++) { - if (matrix.getX(i, j) != matrix.get(j, i)) - return false; - } - - return true; - } - - /** Complex scalar division - real part. */ - private double cdivr; - /** Complex scalar division - imaginary part. */ - private double cdivi; - - /** */ - private void cdiv(double xr, double xi, double yr, double yi) { - double r; - double d; - - if (Math.abs(yr) > Math.abs(yi)) { - r = yi / yr; - d = yr + r * yi; - cdivr = (xr + r * xi) / d; - cdivi = (xi - r * xr) / d; - } - else { - r = yr / yi; - d = yi + r * yr; - cdivr = (r * xr + xi) / d; - cdivi = (r * xi - xr) / d; - } - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/9dccb3db/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/LUDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/LUDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/LUDecomposition.java deleted file mode 100644 index 8b79d9b..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/LUDecomposition.java +++ /dev/null @@ -1,383 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.ignite.ml.math.decompositions; - -import org.apache.ignite.ml.math.Destroyable; -import org.apache.ignite.ml.math.Matrix; -import org.apache.ignite.ml.math.Vector; -import org.apache.ignite.ml.math.exceptions.CardinalityException; -import org.apache.ignite.ml.math.exceptions.SingularMatrixException; - -import static org.apache.ignite.ml.math.util.MatrixUtil.copy; -import static org.apache.ignite.ml.math.util.MatrixUtil.like; -import static org.apache.ignite.ml.math.util.MatrixUtil.likeVector; - -/** - * Calculates the LU-decomposition of a square matrix. - * <p> - * This class is inspired by class from Apache Common Math with similar name.</p> - * - * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html";>MathWorld</a> - * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition";>Wikipedia</a> - * - * <p>TODO: IGNITE-5828, Maybe we should make this class (and other decompositions) Externalizable.</p> - */ -public class LUDecomposition implements Destroyable { - /** Default bound to determine effective singularity in LU decomposition. */ - private static final double DEFAULT_TOO_SMALL = 1e-11; - - /** Pivot permutation associated with LU decomposition. */ - private final Vector pivot; - - /** Parity of the permutation associated with the LU decomposition. */ - private boolean even; - /** Singularity indicator. */ - private boolean singular; - - /** Cached value of L. */ - private Matrix cachedL; - /** Cached value of U. */ - private Matrix cachedU; - /** Cached value of P. */ - private Matrix cachedP; - - /** Original matrix. */ - private Matrix matrix; - - /** Entries of LU decomposition. */ - private Matrix lu; - - /** - * Calculates the LU-decomposition of the given matrix. - * This constructor uses 1e-11 as default value for the singularity - * threshold. - * - * @param matrix Matrix to decompose. - * @throws CardinalityException if matrix is not square. - */ - public LUDecomposition(Matrix matrix) { - this(matrix, DEFAULT_TOO_SMALL); - } - - /** - * Calculates the LUP-decomposition of the given matrix. - * - * @param matrix Matrix to decompose. - * @param singularityThreshold threshold (based on partial row norm). - * @throws CardinalityException if matrix is not square. - */ - public LUDecomposition(Matrix matrix, double singularityThreshold) { - assert matrix != null; - - int rows = matrix.rowSize(); - int cols = matrix.columnSize(); - - if (rows != cols) - throw new CardinalityException(rows, cols); - - this.matrix = matrix; - - lu = copy(matrix); - - pivot = likeVector(matrix); - - for (int i = 0; i < pivot.size(); i++) - pivot.setX(i, i); - - even = true; - singular = false; - - cachedL = null; - cachedU = null; - cachedP = null; - - for (int col = 0; col < cols; col++) { - - //upper - for (int row = 0; row < col; row++) { - Vector luRow = lu.viewRow(row); - double sum = luRow.get(col); - - for (int i = 0; i < row; i++) - sum -= luRow.getX(i) * lu.getX(i, col); - - luRow.setX(col, sum); - } - - // permutation row - int max = col; - - double largest = Double.NEGATIVE_INFINITY; - - // lower - for (int row = col; row < rows; row++) { - Vector luRow = lu.viewRow(row); - double sum = luRow.getX(col); - - for (int i = 0; i < col; i++) - sum -= luRow.getX(i) * lu.getX(i, col); - - luRow.setX(col, sum); - - if (Math.abs(sum) > largest) { - largest = Math.abs(sum); - max = row; - } - } - - // Singularity check - if (Math.abs(lu.getX(max, col)) < singularityThreshold) { - singular = true; - return; - } - - // Pivot if necessary - if (max != col) { - double tmp; - Vector luMax = lu.viewRow(max); - Vector luCol = lu.viewRow(col); - - for (int i = 0; i < cols; i++) { - tmp = luMax.getX(i); - luMax.setX(i, luCol.getX(i)); - luCol.setX(i, tmp); - } - - int temp = (int)pivot.getX(max); - pivot.setX(max, pivot.getX(col)); - pivot.setX(col, temp); - - even = !even; - } - - // Divide the lower elements by the "winning" diagonal elt. - final double luDiag = lu.getX(col, col); - - for (int row = col + 1; row < cols; row++) { - double val = lu.getX(row, col) / luDiag; - lu.setX(row, col, val); - } - } - } - - /** - * Destroys decomposition components and other internal components of decomposition. - */ - @Override public void destroy() { - if (cachedL != null) - cachedL.destroy(); - if (cachedU != null) - cachedU.destroy(); - if (cachedP != null) - cachedP.destroy(); - lu.destroy(); - } - - /** - * Returns the matrix L of the decomposition. - * <p>L is a lower-triangular matrix</p> - * - * @return the L matrix (or null if decomposed matrix is singular). - */ - public Matrix getL() { - if ((cachedL == null) && !singular) { - final int m = pivot.size(); - - cachedL = like(matrix); - cachedL.assign(0.0); - - for (int i = 0; i < m; ++i) { - for (int j = 0; j < i; ++j) - cachedL.setX(i, j, lu.getX(i, j)); - - cachedL.setX(i, i, 1.0); - } - } - - return cachedL; - } - - /** - * Returns the matrix U of the decomposition. - * <p>U is an upper-triangular matrix</p> - * - * @return the U matrix (or null if decomposed matrix is singular). - */ - public Matrix getU() { - if ((cachedU == null) && !singular) { - final int m = pivot.size(); - - cachedU = like(matrix); - cachedU.assign(0.0); - - for (int i = 0; i < m; ++i) - for (int j = i; j < m; ++j) - cachedU.setX(i, j, lu.getX(i, j)); - } - - return cachedU; - } - - /** - * Returns the P rows permutation matrix. - * <p>P is a sparse matrix with exactly one element set to 1.0 in - * each row and each column, all other elements being set to 0.0.</p> - * <p>The positions of the 1 elements are given by the {@link #getPivot() - * pivot permutation vector}.</p> - * - * @return the P rows permutation matrix (or null if decomposed matrix is singular). - * @see #getPivot() - */ - public Matrix getP() { - if ((cachedP == null) && !singular) { - final int m = pivot.size(); - - cachedP = like(matrix); - cachedP.assign(0.0); - - for (int i = 0; i < m; ++i) - cachedP.setX(i, (int)pivot.get(i), 1.0); - } - - return cachedP; - } - - /** - * Returns the pivot permutation vector. - * - * @return the pivot permutation vector. - * @see #getP() - */ - public Vector getPivot() { - return pivot.copy(); - } - - /** - * Return the determinant of the matrix. - * - * @return determinant of the matrix. - */ - public double determinant() { - if (singular) - return 0; - - final int m = pivot.size(); - double determinant = even ? 1 : -1; - - for (int i = 0; i < m; i++) - determinant *= lu.getX(i, i); - - return determinant; - } - - /** - * @param b Vector to solve using this decomposition. - * @return Solution vector. - */ - public Vector solve(Vector b) { - final int m = pivot.size(); - - if (b.size() != m) - throw new CardinalityException(b.size(), m); - - if (singular) - throw new SingularMatrixException(); - - final double[] bp = new double[m]; - - // Apply permutations to b - for (int row = 0; row < m; row++) - bp[row] = b.get((int)pivot.get(row)); - - // Solve LY = b - for (int col = 0; col < m; col++) { - final double bpCol = bp[col]; - - for (int i = col + 1; i < m; i++) - bp[i] -= bpCol * lu.get(i, col); - } - - // Solve UX = Y - for (int col = m - 1; col >= 0; col--) { - bp[col] /= lu.get(col, col); - final double bpCol = bp[col]; - - for (int i = 0; i < col; i++) - bp[i] -= bpCol * lu.get(i, col); - } - - return b.like(m).assign(bp); - } - - /** - * @param b Matrix to solve using this decomposition. - * @return Solution matrix. - */ - public Matrix solve(Matrix b) { - final int m = pivot.size(); - - if (b.rowSize() != m) - throw new CardinalityException(b.rowSize(), m); - - if (singular) - throw new SingularMatrixException(); - - final int nColB = b.columnSize(); - - // Apply permutations to b - final double[][] bp = new double[m][nColB]; - for (int row = 0; row < m; row++) { - final double[] bpRow = bp[row]; - final int pRow = (int)pivot.get(row); - - for (int col = 0; col < nColB; col++) - bpRow[col] = b.get(pRow, col); - } - - // Solve LY = b - for (int col = 0; col < m; col++) { - final double[] bpCol = bp[col]; - for (int i = col + 1; i < m; i++) { - final double[] bpI = bp[i]; - final double luICol = lu.get(i, col); - - for (int j = 0; j < nColB; j++) - bpI[j] -= bpCol[j] * luICol; - } - } - - // Solve UX = Y - for (int col = m - 1; col >= 0; col--) { - final double[] bpCol = bp[col]; - final double luDiag = lu.getX(col, col); - - for (int j = 0; j < nColB; j++) - bpCol[j] /= luDiag; - - for (int i = 0; i < col; i++) { - final double[] bpI = bp[i]; - final double luICol = lu.get(i, col); - - for (int j = 0; j < nColB; j++) - bpI[j] -= bpCol[j] * luICol; - } - } - - return b.like(b.rowSize(), b.columnSize()).assign(bp); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/9dccb3db/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDSolver.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDSolver.java b/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDSolver.java deleted file mode 100644 index bb591ee..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDSolver.java +++ /dev/null @@ -1,197 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.ignite.ml.math.decompositions; - -import java.io.Serializable; -import org.apache.ignite.ml.math.Matrix; -import org.apache.ignite.ml.math.Vector; -import org.apache.ignite.ml.math.exceptions.NoDataException; -import org.apache.ignite.ml.math.exceptions.NullArgumentException; -import org.apache.ignite.ml.math.exceptions.SingularMatrixException; -import org.apache.ignite.ml.math.functions.Functions; -import org.apache.ignite.ml.math.util.MatrixUtil; - -import static org.apache.ignite.ml.math.util.MatrixUtil.like; - -/** - * For an {@code m x n} matrix {@code A} with {@code m >= n}, the QR decomposition - * is an {@code m x n} orthogonal matrix {@code Q} and an {@code n x n} upper - * triangular matrix {@code R} so that {@code A = Q*R}. - */ -public class QRDSolver implements Serializable { - /** */ - private final Matrix q; - - /** */ - private final Matrix r; - - /** - * Constructs a new QR decomposition solver object. - * - * @param q An orthogonal matrix. - * @param r An upper triangular matrix - */ - public QRDSolver(Matrix q, Matrix r) { - this.q = q; - this.r = r; - } - - /** - * Least squares solution of {@code A*X = B}; {@code returns X}. - * - * @param mtx A matrix with as many rows as {@code A} and any number of cols. - * @return {@code X<} that minimizes the two norm of {@code Q*R*X - B}. - * @throws IllegalArgumentException if {@code B.rows() != A.rows()}. - */ - public Matrix solve(Matrix mtx) { - if (mtx.rowSize() != q.rowSize()) - throw new IllegalArgumentException("Matrix row dimensions must agree."); - - int cols = mtx.columnSize(); - Matrix x = like(r, r.columnSize(), cols); - - Matrix qt = q.transpose(); - Matrix y = qt.times(mtx); - - for (int k = Math.min(r.columnSize(), q.rowSize()) - 1; k >= 0; k--) { - // X[k,] = Y[k,] / R[k,k], note that X[k,] starts with 0 so += is same as = - x.viewRow(k).map(y.viewRow(k), Functions.plusMult(1 / r.get(k, k))); - - if (k == 0) - continue; - - // Y[0:(k-1),] -= R[0:(k-1),k] * X[k,] - Vector rCol = r.viewColumn(k).viewPart(0, k); - - for (int c = 0; c < cols; c++) - y.viewColumn(c).viewPart(0, k).map(rCol, Functions.plusMult(-x.get(k, c))); - } - - return x; - } - - /** - * Least squares solution of {@code A*X = B}; {@code returns X}. - * - * @param vec A vector with as many rows as {@code A}. - * @return {@code X<} that minimizes the two norm of {@code Q*R*X - B}. - * @throws IllegalArgumentException if {@code B.rows() != A.rows()}. - */ - public Vector solve(Vector vec) { - if (vec == null) - throw new NullArgumentException(); - if (vec.size() == 0) - throw new NoDataException(); - // TODO: IGNITE-5826, Should we copy here? - - Matrix res = solve(vec.likeMatrix(vec.size(), 1).assignColumn(0, vec)); - - return vec.like(res.rowSize()).assign(res.viewColumn(0)); - } - - /** - * <p>Compute the "hat" matrix. - * </p> - * <p>The hat matrix is defined in terms of the design matrix X - * by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup> - * </p> - * <p>The implementation here uses the QR decomposition to compute the - * hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the - * p-dimensional identity matrix augmented by 0's. This computational - * formula is from "The Hat Matrix in Regression and ANOVA", - * David C. Hoaglin and Roy E. Welsch, - * <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22. - * </p> - * <p>Data for the model must have been successfully loaded using one of - * the {@code newSampleData} methods before invoking this method; otherwise - * a {@code NullPointerException} will be thrown.</p> - * - * @return the hat matrix - * @throws NullPointerException unless method {@code newSampleData} has been called beforehand. - */ - public Matrix calculateHat() { - // Create augmented identity matrix - // No try-catch or advertised NotStrictlyPositiveException - NPE above if n < 3 - Matrix augI = MatrixUtil.like(q, q.columnSize(), q.columnSize()); - - int n = augI.columnSize(); - int p = r.columnSize(); - - for (int i = 0; i < n; i++) - for (int j = 0; j < n; j++) - if (i == j && i < p) - augI.setX(i, j, 1d); - else - augI.setX(i, j, 0d); - - // Compute and return Hat matrix - // No DME advertised - args valid if we get here - return q.times(augI).times(q.transpose()); - } - - /** - * <p>Calculates the variance-covariance matrix of the regression parameters. - * </p> - * <p>Var(b) = (X<sup>T</sup>X)<sup>-1</sup> - * </p> - * <p>Uses QR decomposition to reduce (X<sup>T</sup>X)<sup>-1</sup> - * to (R<sup>T</sup>R)<sup>-1</sup>, with only the top p rows of - * R included, where p = the length of the beta vector.</p> - * - * <p>Data for the model must have been successfully loaded using one of - * the {@code newSampleData} methods before invoking this method; otherwise - * a {@code NullPointerException} will be thrown.</p> - * - * @param p Size of the beta variance-covariance matrix - * @return The beta variance-covariance matrix - * @throws SingularMatrixException if the design matrix is singular - * @throws NullPointerException if the data for the model have not been loaded - */ - public Matrix calculateBetaVariance(int p) { - Matrix rAug = MatrixUtil.copy(r.viewPart(0, p, 0, p)); - Matrix rInv = rAug.inverse(); - - return rInv.times(rInv.transpose()); - } - - /** {@inheritDoc} */ - @Override public boolean equals(Object o) { - if (this == o) - return true; - if (o == null || getClass() != o.getClass()) - return false; - - QRDSolver solver = (QRDSolver)o; - - return q.equals(solver.q) && r.equals(solver.r); - } - - /** {@inheritDoc} */ - @Override public int hashCode() { - int res = q.hashCode(); - res = 31 * res + r.hashCode(); - return res; - } - - /** - * Returns a rough string rendition of a QRD solver. - */ - @Override public String toString() { - return String.format("QRD Solver(%d x %d)", q.rowSize(), r.columnSize()); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/9dccb3db/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDecomposition.java deleted file mode 100644 index c069683..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/ml/math/decompositions/QRDecomposition.java +++ /dev/null @@ -1,212 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.ignite.ml.math.decompositions; - -import org.apache.ignite.ml.math.Destroyable; -import org.apache.ignite.ml.math.Matrix; -import org.apache.ignite.ml.math.Vector; -import org.apache.ignite.ml.math.exceptions.SingularMatrixException; -import org.apache.ignite.ml.math.functions.Functions; - -import static org.apache.ignite.ml.math.util.MatrixUtil.copy; -import static org.apache.ignite.ml.math.util.MatrixUtil.like; - -/** - * For an {@code m x n} matrix {@code A} with {@code m >= n}, the QR decomposition - * is an {@code m x n} orthogonal matrix {@code Q} and an {@code n x n} upper - * triangular matrix {@code R} so that {@code A = Q*R}. - */ -public class QRDecomposition implements Destroyable { - /** */ - private final Matrix q; - /** */ - private final Matrix r; - - /** */ - private final Matrix mType; - /** */ - private final boolean fullRank; - - /** */ - private final int rows; - /** */ - private final int cols; - - /** - * @param v Value to be checked for being an ordinary double. - */ - private void checkDouble(double v) { - if (Double.isInfinite(v) || Double.isNaN(v)) - throw new ArithmeticException("Invalid intermediate result"); - } - - /** - * Constructs a new QR decomposition object computed by Householder reflections. - * Threshold for singularity check used in this case is 0. - * - * @param mtx A rectangular matrix. - */ - public QRDecomposition(Matrix mtx) { - this(mtx, 0.0); - } - - /** - * Constructs a new QR decomposition object computed by Householder reflections. - * - * @param mtx A rectangular matrix. - * @param threshold Value used for detecting singularity of {@code R} matrix in decomposition. - */ - public QRDecomposition(Matrix mtx, double threshold) { - assert mtx != null; - - rows = mtx.rowSize(); - - int min = Math.min(mtx.rowSize(), mtx.columnSize()); - - cols = mtx.columnSize(); - - mType = like(mtx, 1, 1); - - Matrix qTmp = copy(mtx); - - boolean fullRank = true; - - r = like(mtx, min, cols); - - for (int i = 0; i < min; i++) { - Vector qi = qTmp.viewColumn(i); - - double alpha = qi.kNorm(2); - - if (Math.abs(alpha) > Double.MIN_VALUE) - qi.map(Functions.div(alpha)); - else { - checkDouble(alpha); - - fullRank = false; - } - - r.set(i, i, alpha); - - for (int j = i + 1; j < cols; j++) { - Vector qj = qTmp.viewColumn(j); - - double norm = qj.kNorm(2); - - if (Math.abs(norm) > Double.MIN_VALUE) { - double beta = qi.dot(qj); - - r.set(i, j, beta); - - if (j < min) - qj.map(qi, Functions.plusMult(-beta)); - } - else - checkDouble(norm); - } - } - - if (cols > min) - q = qTmp.viewPart(0, rows, 0, min).copy(); - else - q = qTmp; - - verifyNonSingularR(threshold); - - this.fullRank = fullRank; - } - - /** {@inheritDoc} */ - @Override public void destroy() { - q.destroy(); - r.destroy(); - mType.destroy(); - } - - /** - * Gets orthogonal factor {@code Q}. - */ - public Matrix getQ() { - return q; - } - - /** - * Gets triangular factor {@code R}. - */ - public Matrix getR() { - return r; - } - - /** - * Returns whether the matrix {@code A} has full rank. - * - * @return true if {@code R}, and hence {@code A} , has full rank. - */ - public boolean hasFullRank() { - return fullRank; - } - - /** - * Least squares solution of {@code A*X = B}; {@code returns X}. - * - * @param mtx A matrix with as many rows as {@code A} and any number of cols. - * @return {@code X<} that minimizes the two norm of {@code Q*R*X - B}. - * @throws IllegalArgumentException if {@code B.rows() != A.rows()}. - */ - public Matrix solve(Matrix mtx) { - return new QRDSolver(q, r).solve(mtx); - } - - /** - * Least squares solution of {@code A*X = B}; {@code returns X}. - * - * @param vec A vector with as many rows as {@code A}. - * @return {@code X<} that minimizes the two norm of {@code Q*R*X - B}. - * @throws IllegalArgumentException if {@code B.rows() != A.rows()}. - */ - public Vector solve(Vector vec) { - return new QRDSolver(q, r).solve(vec); - } - - /** - * Returns a rough string rendition of a QR. - */ - @Override public String toString() { - return String.format("QR(%d x %d, fullRank=%s)", rows, cols, hasFullRank()); - } - - /** - * Check singularity. - * - * @param min Singularity threshold. - * @throws SingularMatrixException if the matrix is singular and {@code raise} is {@code true}. - */ - private void verifyNonSingularR(double min) { - // TODO: IGNITE-5828, Not a very fast approach for distributed matrices. would be nice if we could independently - // check parts on different nodes for singularity and do fold with 'or'. - - final int len = r.columnSize() > r.rowSize() ? r.rowSize() : r.columnSize(); - for (int i = 0; i < len; i++) { - final double d = r.getX(i, i); - if (Math.abs(d) <= min) - throw new SingularMatrixException("Number is too small (%f, while " + - "threshold is %f). Index of diagonal element is (%d, %d)", d, min, i, i); - - } - } -}
