http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java deleted file mode 100644 index 66fe13c..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java +++ /dev/null @@ -1,923 +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.math.decompositions; - -import org.apache.ignite.math.Matrix; -import org.apache.ignite.math.Vector; -import org.apache.ignite.math.functions.Functions; - -/** - * 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 extends DecompositionSupport { - /** 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; - - /** */ - public EigenDecomposition(Matrix matrix) { - this(matrix, isSymmetric(matrix)); - } - - /** */ - 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 - */ - 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.0); - - // 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/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/LUDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/LUDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/LUDecomposition.java deleted file mode 100644 index 82c90ec..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/LUDecomposition.java +++ /dev/null @@ -1,366 +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.math.decompositions; - -import org.apache.ignite.math.Matrix; -import org.apache.ignite.math.Vector; -import org.apache.ignite.math.exceptions.CardinalityException; -import org.apache.ignite.math.exceptions.SingularMatrixException; - -/** - * Calculates the LU-decomposition of a square matrix. - * - * This class inspired by class from Apache Common Math with similar name. - * - * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html";>MathWorld</a> - * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition";>Wikipedia</a> - */ -public class LUDecomposition extends DecompositionSupport { - /** 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; - } - - /** */ - 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); - } - - /** */ - 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/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/QRDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/QRDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/QRDecomposition.java deleted file mode 100644 index 9608ed5..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/QRDecomposition.java +++ /dev/null @@ -1,186 +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.math.decompositions; - -import org.apache.ignite.math.Matrix; -import org.apache.ignite.math.Vector; -import org.apache.ignite.math.functions.Functions; - -/** - * 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 extends DecompositionSupport { - /** */ - 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. - * - * @param mtx A rectangular matrix. - */ - public QRDecomposition(Matrix mtx) { - 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; - - 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) { - if (mtx.rowSize() != rows) - throw new IllegalArgumentException("Matrix row dimensions must agree."); - - int cols = mtx.columnSize(); - - Matrix x = like(mType, this.cols, cols); - - Matrix qt = getQ().transpose(); - Matrix y = qt.times(mtx); - - Matrix r = getR(); - - for (int k = Math.min(this.cols, rows) - 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))); - - // 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; - } - - /** - * Returns a rough string rendition of a QR. - */ - @Override public String toString() { - return String.format("QR(%d x %d, fullRank=%s)", rows, cols, hasFullRank()); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/SingularValueDecomposition.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/SingularValueDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/SingularValueDecomposition.java deleted file mode 100644 index 75eb206..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/SingularValueDecomposition.java +++ /dev/null @@ -1,620 +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.math.decompositions; - -import org.apache.ignite.math.Algebra; -import org.apache.ignite.math.Matrix; - -/** - * Compute a singular value decomposition (SVD) of {@code (l x k)} matrix {@code m}. - * <p>This decomposition can be thought - * as an extension of {@link EigenDecomposition} to rectangular matrices. The factorization we get is following:</p> - * <p>{@code m = u * s * v^{*}}, where</p> - * <ul><li>{@code u} is a real or complex unitary matrix.</li> - * <li>{@code s} is a rectangular diagonal matrix with non-negative real numbers on diagonal - * (these numbers are singular values of {@code m}).</li> - * <li>{@code v} is a real or complex unitary matrix.</li></ul> - * <p>If {@code m} is real then {@code u} and {@code v} are also real.</p> - * <p>See also: <a href="https://en.wikipedia.org/wiki/Singular_value_decomposition";>Wikipedia article on SVD</a>.</p> - * <p>Note: complex case is currently not supported.</p> - */ -public class SingularValueDecomposition extends DecompositionSupport { - // U and V. - /** */ - private final double[][] u; - /** */ - private final double[][] v; - - /** Singular values. */ - private final double[] s; - - /** Row dimension. */ - private final int m; - /** Column dimension. */ - private final int n; - - /** */ - private Matrix arg; - - /** */ - private boolean transpositionNeeded; - - /** - * Singular value decomposition object. - * - * @param arg A rectangular matrix. - */ - public SingularValueDecomposition(Matrix arg) { - assert arg != null; - - this.arg = arg; - - if (arg.rowSize() < arg.columnSize()) - transpositionNeeded = true; - - double[][] a; - - if (transpositionNeeded) { - // Use the transpose matrix. - m = arg.columnSize(); - n = arg.rowSize(); - - a = new double[m][n]; - - for (int i = 0; i < m; i++) - for (int j = 0; j < n; j++) - a[i][j] = arg.get(j, i); - } - else { - m = arg.rowSize(); - n = arg.columnSize(); - - a = new double[m][n]; - - for (int i = 0; i < m; i++) - for (int j = 0; j < n; j++) - a[i][j] = arg.get(i, j); - } - - int nu = Math.min(m, n); - - s = new double[Math.min(m + 1, n)]; - u = new double[m][nu]; - v = new double[n][n]; - - double[] e = new double[n]; - double[] work = new double[m]; - - int nct = Math.min(m - 1, n); - int nrt = Math.max(0, Math.min(n - 2, m)); - - for (int k = 0; k < Math.max(nct, nrt); k++) { - if (k < nct) { - // Compute the transformation for the k-th column and - // place the k-th diagonal in s[k]. Compute 2-norm of k-th - // column without under/overflow. - s[k] = 0; - - for (int i = k; i < m; i++) - s[k] = Algebra.hypot(s[k], a[i][k]); - - if (s[k] != 0.0) { - if (a[k][k] < 0.0) - s[k] = -s[k]; - - for (int i = k; i < m; i++) - a[i][k] /= s[k]; - - a[k][k] += 1.0; - } - - s[k] = -s[k]; - } - - for (int j = k + 1; j < n; j++) { - if (k < nct && s[k] != 0.0) { - // Apply the transformation. - double t = 0; - - for (int i = k; i < m; i++) - t += a[i][k] * a[i][j]; - - t = -t / a[k][k]; - - for (int i = k; i < m; i++) - a[i][j] += t * a[i][k]; - } - - // Place the k-th row of A into e for the - // subsequent calculation of the row transformation. - e[j] = a[k][j]; - } - - if (k < nct) - // Place the transformation in U for subsequent back - // multiplication. - for (int i = k; i < m; i++) - u[i][k] = a[i][k]; - - if (k < nrt) { - // Compute the k-th row transformation and place the - // k-th super-diagonal in e[k]. - // Compute 2-norm without under/overflow. - e[k] = 0; - - for (int i = k + 1; i < n; i++) - e[k] = Algebra.hypot(e[k], e[i]); - - if (e[k] != 0.0) { - if (e[k + 1] < 0.0) - e[k] = -e[k]; - - for (int i = k + 1; i < n; i++) - e[i] /= e[k]; - - e[k + 1] += 1.0; - } - - e[k] = -e[k]; - - if (k + 1 < m && e[k] != 0.0) { - // Apply the transformation. - for (int i = k + 1; i < m; i++) - work[i] = 0.0; - - for (int j = k + 1; j < n; j++) - for (int i = k + 1; i < m; i++) - work[i] += e[j] * a[i][j]; - - for (int j = k + 1; j < n; j++) { - double t = -e[j] / e[k + 1]; - - for (int i = k + 1; i < m; i++) - a[i][j] += t * work[i]; - } - } - - // Place the transformation in V for subsequent - // back multiplication. - for (int i = k + 1; i < n; i++) - v[i][k] = e[i]; - } - } - - // Set up the final bi-diagonal matrix or order p. - int p = Math.min(n, m + 1); - - if (nct < n) - s[nct] = a[nct][nct]; - - if (m < p) - s[p - 1] = 0.0; - - if (nrt + 1 < p) - e[nrt] = a[nrt][p - 1]; - - e[p - 1] = 0.0; - - // Generate U. - for (int j = nct; j < nu; j++) { - for (int i = 0; i < m; i++) - u[i][j] = 0.0; - - u[j][j] = 1.0; - } - - for (int k = nct - 1; k >= 0; k--) { - if (s[k] != 0.0) { - for (int j = k + 1; j < nu; j++) { - double t = 0; - - for (int i = k; i < m; i++) - t += u[i][k] * u[i][j]; - - t = -t / u[k][k]; - - for (int i = k; i < m; i++) - u[i][j] += t * u[i][k]; - } - - for (int i = k; i < m; i++) - u[i][k] = -u[i][k]; - - u[k][k] = 1.0 + u[k][k]; - - for (int i = 0; i < k - 1; i++) - u[i][k] = 0.0; - } - else { - for (int i = 0; i < m; i++) - u[i][k] = 0.0; - - u[k][k] = 1.0; - } - } - - // Generate V. - for (int k = n - 1; k >= 0; k--) { - if (k < nrt && e[k] != 0.0) { - for (int j = k + 1; j < nu; j++) { - double t = 0; - - for (int i = k + 1; i < n; i++) - t += v[i][k] * v[i][j]; - - t = -t / v[k + 1][k]; - - for (int i = k + 1; i < n; i++) - v[i][j] += t * v[i][k]; - } - } - - for (int i = 0; i < n; i++) - v[i][k] = 0.0; - - v[k][k] = 1.0; - } - - // Main iteration loop for the singular values. - int pp = p - 1; - int iter = 0; - - double eps = Math.pow(2.0, -52.0); - double tiny = Math.pow(2.0, -966.0); - - while (p > 0) { - int k; - - for (k = p - 2; k >= -1; k--) { - if (k == -1) - break; - - if (Math.abs(e[k]) <= tiny + eps * (Math.abs(s[k]) + Math.abs(s[k + 1]))) { - e[k] = 0.0; - - break; - } - } - - int kase; - - if (k == p - 2) - kase = 4; - else { - int ks; - - for (ks = p - 1; ks >= k; ks--) { - if (ks == k) - break; - - double t = - (ks != p ? Math.abs(e[ks]) : 0.) + - (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.); - - if (Math.abs(s[ks]) <= tiny + eps * t) { - s[ks] = 0.0; - - break; - } - } - - if (ks == k) - kase = 3; - else if (ks == p - 1) - kase = 1; - else { - kase = 2; - - k = ks; - } - } - - k++; - - // Perform the task indicated by kase. - switch (kase) { - // Deflate negligible s(p). - case 1: { - double f = e[p - 2]; - - e[p - 2] = 0.0; - - for (int j = p - 2; j >= k; j--) { - double t = Algebra.hypot(s[j], f); - double cs = s[j] / t; - double sn = f / t; - - s[j] = t; - - if (j != k) { - f = -sn * e[j - 1]; - e[j - 1] = cs * e[j - 1]; - } - - for (int i = 0; i < n; i++) { - t = cs * v[i][j] + sn * v[i][p - 1]; - - v[i][p - 1] = -sn * v[i][j] + cs * v[i][p - 1]; - v[i][j] = t; - } - } - } - - break; - - // Split at negligible s(k). - case 2: { - double f = e[k - 1]; - e[k - 1] = 0.0; - - for (int j = k; j < p; j++) { - double t = Algebra.hypot(s[j], f); - double cs = s[j] / t; - double sn = f / t; - - s[j] = t; - f = -sn * e[j]; - e[j] = cs * e[j]; - - for (int i = 0; i < m; i++) { - t = cs * u[i][j] + sn * u[i][k - 1]; - - u[i][k - 1] = -sn * u[i][j] + cs * u[i][k - 1]; - u[i][j] = t; - } - } - } - - break; - - // Perform one qr step. - case 3: { - // Calculate the shift. - double scale = Math.max(Math.max(Math.max(Math.max( - Math.abs(s[p - 1]), Math.abs(s[p - 2])), Math.abs(e[p - 2])), - Math.abs(s[k])), Math.abs(e[k])); - - double sp = s[p - 1] / scale; - double spm1 = s[p - 2] / scale; - double epm1 = e[p - 2] / scale; - double sk = s[k] / scale; - double ek = e[k] / scale; - double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0; - double c = sp * epm1 * sp * epm1; - double shift = 0.0; - - if (b != 0.0 || c != 0.0) { - shift = Math.sqrt(b * b + c); - - if (b < 0.0) - shift = -shift; - - shift = c / (b + shift); - } - - double f = (sk + sp) * (sk - sp) + shift; - double g = sk * ek; - - // Chase zeros. - for (int j = k; j < p - 1; j++) { - double t = Algebra.hypot(f, g); - double cs = f / t; - double sn = g / t; - - if (j != k) - e[j - 1] = t; - - f = cs * s[j] + sn * e[j]; - e[j] = cs * e[j] - sn * s[j]; - g = sn * s[j + 1]; - s[j + 1] = cs * s[j + 1]; - - for (int i = 0; i < n; i++) { - t = cs * v[i][j] + sn * v[i][j + 1]; - - v[i][j + 1] = -sn * v[i][j] + cs * v[i][j + 1]; - v[i][j] = t; - } - - t = Algebra.hypot(f, g); - cs = f / t; - sn = g / t; - s[j] = t; - f = cs * e[j] + sn * s[j + 1]; - s[j + 1] = -sn * e[j] + cs * s[j + 1]; - g = sn * e[j + 1]; - e[j + 1] = cs * e[j + 1]; - - if (j < m - 1) - for (int i = 0; i < m; i++) { - t = cs * u[i][j] + sn * u[i][j + 1]; - - u[i][j + 1] = -sn * u[i][j] + cs * u[i][j + 1]; - u[i][j] = t; - } - } - - e[p - 2] = f; - iter = iter + 1; - } - - break; - - // Convergence. - case 4: { - // Make the singular values positive. - if (s[k] <= 0.0) { - s[k] = s[k] < 0.0 ? -s[k] : 0.0; - - for (int i = 0; i <= pp; i++) - v[i][k] = -v[i][k]; - } - - // Order the singular values. - while (k < pp) { - if (s[k] >= s[k + 1]) - break; - - double t = s[k]; - - s[k] = s[k + 1]; - s[k + 1] = t; - - if (k < n - 1) - for (int i = 0; i < n; i++) { - t = v[i][k + 1]; - - v[i][k + 1] = v[i][k]; - v[i][k] = t; - } - - if (k < m - 1) - for (int i = 0; i < m; i++) { - t = u[i][k + 1]; - - u[i][k + 1] = u[i][k]; - u[i][k] = t; - } - - k++; - } - - iter = 0; - p--; - } - - break; - - default: - throw new IllegalStateException(); - } - } - } - - /** - * Gets the two norm condition number, which is {@code max(S) / min(S)} . - */ - public double cond() { - return s[0] / s[Math.min(m, n) - 1]; - } - - /** - * @return the diagonal matrix of singular values. - */ - public Matrix getS() { - double[][] s = new double[n][n]; - - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) - s[i][j] = 0.0; - - s[i][i] = this.s[i]; - } - - return like(arg, n, n).assign(s); - } - - /** - * Gets the diagonal of {@code S}, which is a one-dimensional array of - * singular values. - * - * @return diagonal of {@code S}. - */ - public double[] getSingularValues() { - return s; - } - - /** - * Gets the left singular vectors {@code U}. - * - * @return {@code U} - */ - public Matrix getU() { - if (transpositionNeeded) - return like(arg, v.length, v.length).assign(v); - else { - int numCols = Math.min(m + 1, n); - - Matrix r = like(arg, m, numCols); - - for (int i = 0; i < m; i++) - for (int j = 0; j < numCols; j++) - r.set(i, j, u[i][j]); - - return r; - } - } - - /** - * Gets the right singular vectors {@code V}. - * - * @return {@code V} - */ - public Matrix getV() { - if (transpositionNeeded) { - int numCols = Math.min(m + 1, n); - - Matrix r = like(arg, m, numCols); - - for (int i = 0; i < m; i++) - for (int j = 0; j < numCols; j++) - r.set(i, j, u[i][j]); - - return r; - } - else - return like(arg, v.length, v.length).assign(v); - } - - /** - * Gets the two norm, which is {@code max(S)}. - */ - public double norm2() { - return s[0]; - } - - /** - * Gets effective numerical matrix rank. - */ - public int rank() { - double eps = Math.pow(2.0, -52.0); - double tol = Math.max(m, n) * s[0] * eps; - int r = 0; - - for (double value : s) - if (value > tol) - r++; - - return r; - } - - /** - * Gets [n à n] covariance matrix. - * - * @param minSingularVal Value below which singular values are ignored. - */ - Matrix getCovariance(double minSingularVal) { - Matrix j = like(arg, s.length, s.length); - Matrix vMat = like(arg, v.length, v.length).assign(v); - - for (int i = 0; i < s.length; i++) - j.set(i, i, s[i] >= minSingularVal ? 1 / (s[i] * s[i]) : 0.0); - - return vMat.times(j).times(vMat.transpose()); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/package-info.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/package-info.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/package-info.java deleted file mode 100644 index dcfa0f8..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/package-info.java +++ /dev/null @@ -1,22 +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 description. --> - * Contains matrix decompositions for distributed code algebra. - */ -package org.apache.ignite.math.decompositions; http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/CardinalityException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/CardinalityException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/CardinalityException.java deleted file mode 100644 index fc87a27..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/CardinalityException.java +++ /dev/null @@ -1,38 +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.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * Indicates a cardinality mismatch in matrix or vector operations. - */ -public class CardinalityException extends IgniteException { - /** */ - private static final long serialVersionUID = 0L; - - /** - * Creates new cardinality violation exception. - * - * @param exp Expected cardinality. - * @param act Actual cardinality. - */ - public CardinalityException(int exp, int act) { - super("Cardinality violation [expected=" + exp + ", actual=" + act + "]"); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/ColumnIndexException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/ColumnIndexException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/ColumnIndexException.java deleted file mode 100644 index 7670caf..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/ColumnIndexException.java +++ /dev/null @@ -1,35 +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.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * This exception is used to indicate any error condition accessing matrix elements by invalid column index. - */ -public class ColumnIndexException extends IgniteException { - /** */ - private static final long serialVersionUID = 0L; - - /** - * @param idx Index value that caused this exception. - */ - public ColumnIndexException(int idx) { - super("Invalid (out of bound) column index: " + idx); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/IndexException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/IndexException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/IndexException.java deleted file mode 100644 index 9ada706..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/IndexException.java +++ /dev/null @@ -1,35 +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.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * Indicates an invalid, i.e. out of bound, index on matrix or vector operations. - */ -public class IndexException extends IgniteException { - /** */ - private static final long serialVersionUID = 0L; - - /** - * @param idx Index value that caused this exception. - */ - public IndexException(int idx) { - super("Invalid (out of bound) index: " + idx); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonPositiveDefiniteMatrixException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonPositiveDefiniteMatrixException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonPositiveDefiniteMatrixException.java deleted file mode 100644 index b6017c2..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonPositiveDefiniteMatrixException.java +++ /dev/null @@ -1,20 +0,0 @@ -package org.apache.ignite.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * This exception is used to indicate error condition of matrix elements failing the positivity check. - */ -public class NonPositiveDefiniteMatrixException extends IgniteException { - /** - * Construct an exception. - * - * @param wrong Value that fails the positivity check. - * @param idx Row (and column) index. - * @param threshold Absolute positivity threshold. - */ - public NonPositiveDefiniteMatrixException(double wrong, int idx, double threshold) { - super("Matrix must be positive, wrong element located on diagonal with index " - + idx + " and has value " + wrong + " with this threshold " + threshold); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonSymmetricMatrixException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonSymmetricMatrixException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonSymmetricMatrixException.java deleted file mode 100644 index 8b4cbdb..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/NonSymmetricMatrixException.java +++ /dev/null @@ -1,18 +0,0 @@ -package org.apache.ignite.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * This exception is used to indicate error condition of matrix failing the symmetry check. - */ -public class NonSymmetricMatrixException extends IgniteException { - /** - * @param row Row. - * @param col Column. - * @param threshold Threshold. - */ - public NonSymmetricMatrixException(int row, int col, double threshold) { - super("Symmetric matrix expected, the symmetry is broken on row " - + row + " and col " + col + " with this threshold " + threshold); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/RowIndexException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/RowIndexException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/RowIndexException.java deleted file mode 100644 index f74ae2c..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/RowIndexException.java +++ /dev/null @@ -1,35 +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.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * This exception is used to indicate any error condition accessing matrix elements by invalid row index. - */ -public class RowIndexException extends IgniteException { - /** */ - private static final long serialVersionUID = 0L; - - /** - * @param idx Index value that caused this exception. - */ - public RowIndexException(int idx) { - super("Invalid (out of bound) row index: " + idx); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/SingularMatrixException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/SingularMatrixException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/SingularMatrixException.java deleted file mode 100644 index 4ed3410..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/SingularMatrixException.java +++ /dev/null @@ -1,30 +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.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * Exception to be thrown when a non-singular matrix is expected. - */ -public class SingularMatrixException extends IgniteException { - /** */ - public SingularMatrixException() { - super("Regular (or non-singular) matrix expected."); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnknownProviderException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnknownProviderException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnknownProviderException.java deleted file mode 100644 index 3e6498a..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnknownProviderException.java +++ /dev/null @@ -1,35 +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.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * Indicates that no provider has been found for a given vector or matrix flavor. - */ -public class UnknownProviderException extends IgniteException { - /** */ - private static final long serialVersionUID = 0L; - - /** - * @param flv Flavor (a.k.a. operation performance hints) that has no registered provider for. - */ - public UnknownProviderException(String flv) { - super("No provider has been found for the flavor: " + flv); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnsupportedOperationException.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnsupportedOperationException.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnsupportedOperationException.java deleted file mode 100644 index be5264c..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/UnsupportedOperationException.java +++ /dev/null @@ -1,44 +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.math.exceptions; - -import org.apache.ignite.IgniteException; - -/** - * Indicate that a specific operation is not supported by the underlying implementation. - * In some cases, an operation may be unsupported only in certain cases where, for example, - * it could not be deterministically completed in polynomial time. - */ -public class UnsupportedOperationException extends IgniteException { - /** */ - private static final long serialVersionUID = 0L; - - /** - * @param errMsg Error message. - */ - public UnsupportedOperationException(String errMsg) { - super(errMsg); - } - - /** - * - */ - public UnsupportedOperationException() { - this("Unsupported operation."); - } -} http://git-wip-us.apache.org/repos/asf/ignite/blob/d78e071a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/package-info.java ---------------------------------------------------------------------- diff --git a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/package-info.java b/modules/ml/src/main/java/org/apache/ignite/math/exceptions/package-info.java deleted file mode 100644 index 83f3fa4..0000000 --- a/modules/ml/src/main/java/org/apache/ignite/math/exceptions/package-info.java +++ /dev/null @@ -1,22 +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 description. --> - * Contains exceptions for distributed code algebra. - */ -package org.apache.ignite.math.exceptions; \ No newline at end of file
