Support implicit-Krylov-approximation-based efficient SST
Project: http://git-wip-us.apache.org/repos/asf/incubator-hivemall/repo Commit: http://git-wip-us.apache.org/repos/asf/incubator-hivemall/commit/998203d5 Tree: http://git-wip-us.apache.org/repos/asf/incubator-hivemall/tree/998203d5 Diff: http://git-wip-us.apache.org/repos/asf/incubator-hivemall/diff/998203d5 Branch: refs/heads/JIRA-22/pr-356 Commit: 998203d5e8623d6282c2b187df24e4da7d41c16b Parents: 2bfd127 Author: Takuya Kitazawa <[email protected]> Authored: Wed Sep 28 19:49:48 2016 +0900 Committer: Takuya Kitazawa <[email protected]> Committed: Wed Sep 28 19:49:48 2016 +0900 ---------------------------------------------------------------------- .../anomaly/SingularSpectrumTransform.java | 103 ++++++++-- .../anomaly/SingularSpectrumTransformUDF.java | 27 +++ .../java/hivemall/utils/math/MatrixUtils.java | 203 +++++++++++++++++++ .../anomaly/SingularSpectrumTransformTest.java | 61 ++++-- .../hivemall/utils/math/MatrixUtilsTest.java | 67 ++++++ 5 files changed, 434 insertions(+), 27 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/incubator-hivemall/blob/998203d5/core/src/main/java/hivemall/anomaly/SingularSpectrumTransform.java ---------------------------------------------------------------------- diff --git a/core/src/main/java/hivemall/anomaly/SingularSpectrumTransform.java b/core/src/main/java/hivemall/anomaly/SingularSpectrumTransform.java index c964129..f9f6222 100644 --- a/core/src/main/java/hivemall/anomaly/SingularSpectrumTransform.java +++ b/core/src/main/java/hivemall/anomaly/SingularSpectrumTransform.java @@ -18,9 +18,11 @@ package hivemall.anomaly; import hivemall.anomaly.SingularSpectrumTransformUDF.SingularSpectrumTransformInterface; +import hivemall.anomaly.SingularSpectrumTransformUDF.ScoreFunction; import hivemall.anomaly.SingularSpectrumTransformUDF.Parameters; import hivemall.utils.collections.DoubleRingBuffer; -import org.apache.commons.math3.linear.MatrixUtils; +import hivemall.utils.math.MatrixUtils; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.SingularValueDecomposition; import org.apache.hadoop.hive.ql.metadata.HiveException; @@ -28,6 +30,8 @@ import org.apache.hadoop.hive.serde2.objectinspector.PrimitiveObjectInspector; import org.apache.hadoop.hive.serde2.objectinspector.primitive.PrimitiveObjectInspectorUtils; import java.util.Arrays; +import java.util.TreeMap; +import java.util.Collections; import javax.annotation.Nonnull; @@ -37,6 +41,9 @@ final class SingularSpectrumTransform implements SingularSpectrumTransformInterf private final PrimitiveObjectInspector oi; @Nonnull + private final ScoreFunction scoreFunc; + + @Nonnull private final int window; @Nonnull private final int nPastWindow; @@ -50,15 +57,22 @@ final class SingularSpectrumTransform implements SingularSpectrumTransformInterf private final int currentOffset; @Nonnull private final int r; + @Nonnull + private final int k; @Nonnull private final DoubleRingBuffer xRing; @Nonnull private final double[] xSeries; + @Nonnull + private final double[] q; + SingularSpectrumTransform(@Nonnull Parameters params, @Nonnull PrimitiveObjectInspector oi) { this.oi = oi; + this.scoreFunc = params.scoreFunc; + this.window = params.w; this.nPastWindow = params.n; this.nCurrentWindow = params.m; @@ -66,6 +80,7 @@ final class SingularSpectrumTransform implements SingularSpectrumTransformInterf this.currentSize = window + nCurrentWindow; this.currentOffset = params.g; this.r = params.r; + this.k = params.k; // (w + n) past samples for the n-past-windows // (w + m) current samples for the m-current-windows, starting from offset g @@ -74,6 +89,18 @@ final class SingularSpectrumTransform implements SingularSpectrumTransformInterf this.xRing = new DoubleRingBuffer(holdSampleSize); this.xSeries = new double[holdSampleSize]; + + this.q = new double[window]; + double norm = 0.d; + for (int i = 0; i < window; i++) { + this.q[i] = Math.random(); + norm += q[i] * q[i]; + } + norm = Math.sqrt(norm); + // normalize + for (int i = 0; i < window; i++) { + this.q[i] = q[i] / norm; + } } @Override @@ -86,25 +113,39 @@ final class SingularSpectrumTransform implements SingularSpectrumTransformInterf if (!xRing.isFull()) { outScores[0] = 0.d; } else { - outScores[0] = computeScore(); + // create past trajectory matrix and find its left singular vectors + RealMatrix H = new Array2DRowRealMatrix(new double[window][nPastWindow]); + for (int i = 0; i < nPastWindow; i++) { + H.setColumn(i, Arrays.copyOfRange(xSeries, i, i + window)); + } + + // create current trajectory matrix and find its left singular vectors + RealMatrix G = new Array2DRowRealMatrix(new double[window][nCurrentWindow]); + int currentHead = pastSize + currentOffset; + for (int i = 0; i < nCurrentWindow; i++) { + G.setColumn(i, Arrays.copyOfRange(xSeries, currentHead + i, currentHead + i + window)); + } + + switch (scoreFunc) { + case svd: + outScores[0] = computeScoreSVD(H, G); + break; + case ika: + outScores[0] = computeScoreIKA(H, G); + break; + default: + throw new IllegalStateException("Unexpected score function: " + scoreFunc); + } } } - private double computeScore() { - // create past trajectory matrix and find its left singular vectors - RealMatrix H = MatrixUtils.createRealMatrix(window, nPastWindow); - for (int i = 0; i < nPastWindow; i++) { - H.setColumn(i, Arrays.copyOfRange(xSeries, i, i + window)); - } + /** + * Singular Value Decomposition (SVD) based naive scoring. + */ + private double computeScoreSVD(@Nonnull final RealMatrix H, @Nonnull final RealMatrix G) { SingularValueDecomposition svdH = new SingularValueDecomposition(H); RealMatrix UT = svdH.getUT(); - // create current trajectory matrix and find its left singular vectors - RealMatrix G = MatrixUtils.createRealMatrix(window, nCurrentWindow); - int currentHead = pastSize + currentOffset; - for (int i = 0; i < nCurrentWindow; i++) { - G.setColumn(i, Arrays.copyOfRange(xSeries, currentHead + i, currentHead + i + window)); - } SingularValueDecomposition svdG = new SingularValueDecomposition(G); RealMatrix Q = svdG.getU(); @@ -115,4 +156,38 @@ final class SingularSpectrumTransform implements SingularSpectrumTransformInterf return 1.d - s[0]; } + + /** + * Implicit Krylov Approximation (IKA) based naive scoring. + * + * Number of iterations for the Power method and QR method is fixed to 1 for efficiency. + * This may cause failure (i.e. meaningless scores) depending on datasets and initial values. + * + */ + private double computeScoreIKA(@Nonnull final RealMatrix H, @Nonnull final RealMatrix G) { + // assuming n = m = window, and keep track the left singular vector as `q` + double firstSingularValue = MatrixUtils.power1(G, q, 1, q, new double[window]); + + RealMatrix T = new Array2DRowRealMatrix(new double[k][k]); + MatrixUtils.lanczosTridiagonalization(H.multiply(H.transpose()), q, T); + + double[] eigvals = new double[k]; + RealMatrix eigvecs = new Array2DRowRealMatrix(new double[k][k]); + MatrixUtils.tridiagonalEigen(T, 1, eigvals, eigvecs); + + // tridiagonalEigen() returns unordered eigenvalues, + // so the top-r eigenvectors should be picked carefully + TreeMap<Double, Integer> map = new TreeMap<Double, Integer>(Collections.reverseOrder()); + for (int i = 0; i < k; i++) { + map.put(eigvals[i], i); + } + Object[] sortedIndices = map.values().toArray(); + + double s = 0.d; + for (int i = 0; i < r; i++) { + double v = eigvecs.getEntry(0, (int)sortedIndices[i]); + s += v * v; + } + return 1.d - Math.sqrt(s); + } } http://git-wip-us.apache.org/repos/asf/incubator-hivemall/blob/998203d5/core/src/main/java/hivemall/anomaly/SingularSpectrumTransformUDF.java ---------------------------------------------------------------------- diff --git a/core/src/main/java/hivemall/anomaly/SingularSpectrumTransformUDF.java b/core/src/main/java/hivemall/anomaly/SingularSpectrumTransformUDF.java index 64b7d20..5f8633d 100644 --- a/core/src/main/java/hivemall/anomaly/SingularSpectrumTransformUDF.java +++ b/core/src/main/java/hivemall/anomaly/SingularSpectrumTransformUDF.java @@ -89,6 +89,8 @@ public final class SingularSpectrumTransformUDF extends UDFWithOptions { "Number of singular vectors (i.e. principal components) [default: 3]"); opts.addOption("k", "n_dim", true, "Number of dimensions for the Krylov subspaces [default: 5 (`2*r` if `r` is even, `2*r-1` otherwise)]"); + opts.addOption("score", "scorefunc", true, + "Score function [default: svd, ika]"); opts.addOption("th", "threshold", true, "Score threshold (inclusive) for determining change-point existence [default: -1, do not output decision]"); return opts; @@ -105,6 +107,12 @@ public final class SingularSpectrumTransformUDF extends UDFWithOptions { this._params.r = Primitives.parseInt(cl.getOptionValue("r"), _params.r); this._params.k = Primitives.parseInt( cl.getOptionValue("k"), (_params.r % 2 == 0) ? (2 * _params.r) : (2 * _params.r - 1)); + + this._params.scoreFunc = ScoreFunction.resolve(cl.getOptionValue("scorefunc", ScoreFunction.svd.name())); + if ((_params.w != _params.n || _params.w != _params.m) && _params.scoreFunc == ScoreFunction.ika) { + throw new UDFArgumentException("IKA-based efficient SST requires w = n = m"); + } + this._params.changepointThreshold = Primitives.parseDouble( cl.getOptionValue("th"), _params.changepointThreshold); @@ -196,13 +204,32 @@ public final class SingularSpectrumTransformUDF extends UDFWithOptions { int g = -30; int r = 3; int k = 5; + ScoreFunction scoreFunc = ScoreFunction.svd; double changepointThreshold = -1.d; Parameters() {} + + void set(@Nonnull ScoreFunction func) { + this.scoreFunc = func; + } } public interface SingularSpectrumTransformInterface { void update(@Nonnull Object arg, @Nonnull double[] outScores) throws HiveException; } + public enum ScoreFunction { + svd, ika; + + static ScoreFunction resolve(@Nullable final String name) { + if (svd.name().equalsIgnoreCase(name)) { + return svd; + } else if (ika.name().equalsIgnoreCase(name)) { + return ika; + } else { + throw new IllegalArgumentException("Unsupported ScoreFunction: " + name); + } + } + } + } http://git-wip-us.apache.org/repos/asf/incubator-hivemall/blob/998203d5/core/src/main/java/hivemall/utils/math/MatrixUtils.java ---------------------------------------------------------------------- diff --git a/core/src/main/java/hivemall/utils/math/MatrixUtils.java b/core/src/main/java/hivemall/utils/math/MatrixUtils.java index 840df41..aaf9d4a 100644 --- a/core/src/main/java/hivemall/utils/math/MatrixUtils.java +++ b/core/src/main/java/hivemall/utils/math/MatrixUtils.java @@ -26,11 +26,16 @@ import org.apache.commons.math3.linear.BlockRealMatrix; import org.apache.commons.math3.linear.DecompositionSolver; import org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor; import org.apache.commons.math3.linear.LUDecomposition; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealMatrixPreservingVisitor; import org.apache.commons.math3.linear.SingularMatrixException; import org.apache.commons.math3.linear.SingularValueDecomposition; +import java.util.Arrays; + public final class MatrixUtils { private MatrixUtils() {} @@ -493,4 +498,202 @@ public final class MatrixUtils { return A; } + /** + * Find the first singular vector/value of a matrix A based on the Power method. + * + * http://www.cs.yale.edu/homes/el327/datamining2013aFiles/07_singular_value_decomposition.pdf + * + * @param A target matrix + * @param x0 initial vector + * @param nIter number of iterations for the Power method + * @param u 1st left singular vector + * @param v 1st right singular vector + * @return 1st singular value + */ + @Nonnull + public static double power1(@Nonnull final RealMatrix A, @Nonnull final double[] x0, final int nIter, + @Nonnull final double[] u, @Nonnull final double[] v) { + Preconditions.checkArgument(A.getColumnDimension() == x0.length, + "Column size of A and length of x should be same"); + Preconditions.checkArgument(A.getRowDimension() == u.length, + "Row size of A and length of u should be same"); + Preconditions.checkArgument(x0.length == v.length, "Length of x and u should be same"); + Preconditions.checkArgument(nIter >= 1, "Invalid number of iterations: " + nIter); + + RealMatrix AtA = A.transpose().multiply(A); + + RealVector x = new ArrayRealVector(x0); + for (int i = 0; i < nIter; i++) { + x = AtA.operate(x); + } + + double xNorm = x.getNorm(); + for (int i = 0, n = v.length; i < n; i++) { + v[i] = x.getEntry(i) / xNorm; + } + + RealVector Av = new ArrayRealVector(A.operate(v)); + double s = Av.getNorm(); + + for (int i = 0, n = u.length; i < n; i++) { + u[i] = Av.getEntry(i) / s; + } + + return s; + } + + /** + * Lanczos tridiagonalization for a symmetric matrix C to make s * s tridiagonal matrix T. + * + * http://www.cas.mcmaster.ca/~qiao/publications/spie05.pdf + * + * @param C target symmetric matrix + * @param a initial vector + * @param T result is stored here + */ + @Nonnull + public static void lanczosTridiagonalization(@Nonnull final RealMatrix C, @Nonnull final double[] a, + @Nonnull final RealMatrix T) { + Preconditions.checkArgument(Arrays.deepEquals(C.getData(), C.transpose().getData()), + "Target matrix C must be a symmetric matrix"); + Preconditions.checkArgument(C.getColumnDimension() == a.length, + "Column size of A and length of a should be same"); + Preconditions.checkArgument(T.getRowDimension() == T.getColumnDimension(), + "T must be a square matrix"); + + int s = T.getRowDimension(); + + // initialize T with zeros + T.setSubMatrix(new double[s][s], 0, 0); + + RealVector a0 = new ArrayRealVector(new double[a.length]); + RealVector r = new ArrayRealVector(a); + + double beta0 = 1.d; + + for (int i = 0; i < s; i++) { + RealVector a1 = r.mapDivide(beta0); + RealVector Ca1 = C.operate(a1); + + double alpha1 = a1.dotProduct(Ca1); + + r = Ca1.add(a1.mapMultiply(-1.d * alpha1)).add(a0.mapMultiply(-1.d * beta0)); + + double beta1 = r.getNorm(); + + T.setEntry(i, i, alpha1); + if (i - 1 >= 0) { + T.setEntry(i, i - 1, beta0); + } + if (i + 1 < s) { + T.setEntry(i, i + 1, beta1); + } + + a0 = a1.copy(); + beta0 = beta1; + } + } + + /** + * QR decomposition for a tridiagonal matrix T. + * + * https://gist.github.com/lightcatcher/8118181 + * http://www.ericmart.in/blog/optimizing_julia_tridiag_qr + * + * @param T target tridiagonal matrix + * @param R output matrix for R which is the same shape as T + * @param Qt output matrix for Q.T which is the same shape an T + */ + @Nonnull + public static void tridiagonalQR(@Nonnull final RealMatrix T, + @Nonnull final RealMatrix R, @Nonnull final RealMatrix Qt) { + int n = T.getRowDimension(); + Preconditions.checkArgument(n == R.getRowDimension() && n == R.getColumnDimension(), + "T and R must be the same shape"); + Preconditions.checkArgument(n == Qt.getRowDimension() && n == Qt.getColumnDimension(), + "T and Qt must be the same shape"); + + // initial R = T + R.setSubMatrix(T.getData(), 0, 0); + + // initial Qt = identity + Qt.setSubMatrix(new double[n][n], 0, 0); + for (int i = 0; i < n; i++) { + Qt.setEntry(i, i, 1); + } + + for (int i = 0; i < n - 1; i++) { + // Householder projection for a vector x + // https://en.wikipedia.org/wiki/Householder_transformation + RealVector x = T.getSubMatrix(i, i + 1, i, i).getColumnVector(0); + + double x0 = x.getEntry(0); + double sign = 0.d; + if (x0 < 0.d) { + sign = -1.d; + } else if (x0 > 0.d) { + sign = 1.d; + } + + x.setEntry(0, x0 + sign * x.getNorm()); + x = x.unitVector(); + + RealMatrix subR = R.getSubMatrix(i, i + 1, 0, n - 1); + R.setSubMatrix(subR.subtract(x.outerProduct(subR.preMultiply(x)).scalarMultiply(2)).getData(), i, 0); + + RealMatrix subQt = Qt.getSubMatrix(i, i + 1, 0, n - 1); + Qt.setSubMatrix(subQt.subtract(x.outerProduct(subQt.preMultiply(x)).scalarMultiply(2)).getData(), i, 0); + } + } + + /** + * Find eigenvalues and eigenvectors of given tridiagonal matrix T. + * + * http://web.csulb.edu/~tgao/math423/s94.pdf + * http://stats.stackexchange.com/questions/20643/finding-matrix-eigenvectors-using-qr-decomposition + * + * @param T target tridiagonal matrix + * @param nIter number of iterations for the QR method + * @param eigvals eigenvalues are stored here + * @param eigvecs eigenvectors are stored here + */ + @Nonnull + public static void tridiagonalEigen(@Nonnull final RealMatrix T, @Nonnull final int nIter, + @Nonnull final double[] eigvals, @Nonnull final RealMatrix eigvecs) { + Preconditions.checkArgument(Arrays.deepEquals(T.getData(), T.transpose().getData()), + "Target matrix T must be a symmetric (tridiagonal) matrix"); + Preconditions.checkArgument(eigvecs.getRowDimension() == eigvecs.getColumnDimension(), + "eigvecs must be a square matrix"); + Preconditions.checkArgument(T.getRowDimension() == eigvecs.getRowDimension(), + "T and eigvecs must be the same shape"); + Preconditions.checkArgument(eigvals.length == eigvecs.getRowDimension(), + "Number of eigenvalues and eigenvectors must be same"); + + int nEig = eigvals.length; + + // initialize eigvecs as an identity matrix + eigvecs.setSubMatrix(new double[nEig][nEig], 0, 0); + for (int i = 0; i < nEig; i++) { + eigvecs.setEntry(i, i, 1); + } + + RealMatrix T_ = T.copy(); + + for (int i = 0; i < nIter; i++) { + // QR decomposition for the tridiagonal matrix T + RealMatrix R = new Array2DRowRealMatrix(new double[nEig][nEig]); + RealMatrix Qt = new Array2DRowRealMatrix(new double[nEig][nEig]); + tridiagonalQR(T_, R, Qt); + + RealMatrix Q = Qt.transpose(); + T_ = R.multiply(Q); + eigvecs.setSubMatrix(eigvecs.multiply(Q).getData(), 0, 0); + } + + // diagonal elements correspond to the eigenvalues + for (int i = 0; i < nEig; i++) { + eigvals[i] = T_.getEntry(i, i); + } + } + } http://git-wip-us.apache.org/repos/asf/incubator-hivemall/blob/998203d5/core/src/test/java/hivemall/anomaly/SingularSpectrumTransformTest.java ---------------------------------------------------------------------- diff --git a/core/src/test/java/hivemall/anomaly/SingularSpectrumTransformTest.java b/core/src/test/java/hivemall/anomaly/SingularSpectrumTransformTest.java index d4f119f..44d114d 100644 --- a/core/src/test/java/hivemall/anomaly/SingularSpectrumTransformTest.java +++ b/core/src/test/java/hivemall/anomaly/SingularSpectrumTransformTest.java @@ -17,6 +17,7 @@ */ package hivemall.anomaly; +import hivemall.anomaly.SingularSpectrumTransformUDF.ScoreFunction; import hivemall.anomaly.SingularSpectrumTransformUDF.Parameters; import java.io.BufferedReader; @@ -37,8 +38,45 @@ public class SingularSpectrumTransformTest { private static final boolean DEBUG = false; @Test - public void testSST() throws IOException, HiveException { + public void testSVDSST() throws IOException, HiveException { + int numChangepoints = detectSST(ScoreFunction.svd, 0.95d); + Assert.assertTrue("#changepoints SHOULD be greater than 0: " + numChangepoints, + numChangepoints > 0); + Assert.assertTrue("#changepoints SHOULD be less than 5: " + numChangepoints, + numChangepoints < 5); + } + + @Test + public void testIKASST() throws IOException, HiveException { + int numChangepoints = detectSST(ScoreFunction.ika, 0.65d); + Assert.assertTrue("#changepoints SHOULD be greater than 0: " + numChangepoints, + numChangepoints > 0); + Assert.assertTrue("#changepoints SHOULD be less than 5: " + numChangepoints, + numChangepoints < 5); + } + + @Test + public void testSVDTwitterData() throws IOException, HiveException { + int numChangepoints = detectTwitterData(ScoreFunction.svd, 0.005d); + Assert.assertTrue("#changepoints SHOULD be greater than 0: " + numChangepoints, + numChangepoints > 0); + Assert.assertTrue("#changepoints SHOULD be less than 5: " + numChangepoints, + numChangepoints < 5); + } + + @Test + public void testIKATwitterData() throws IOException, HiveException { + int numChangepoints = detectTwitterData(ScoreFunction.ika, 0.0175d); + Assert.assertTrue("#changepoints SHOULD be greater than 0: " + numChangepoints, + numChangepoints > 0); + Assert.assertTrue("#changepoints SHOULD be less than 5: " + numChangepoints, + numChangepoints < 5); + } + + private static int detectSST(@Nonnull final ScoreFunction scoreFunc, + @Nonnull final double threshold) throws IOException, HiveException { Parameters params = new Parameters(); + params.set(scoreFunc); PrimitiveObjectInspector oi = PrimitiveObjectInspectorFactory.javaDoubleObjectInspector; SingularSpectrumTransform sst = new SingularSpectrumTransform(params, oi); double[] outScores = new double[1]; @@ -51,19 +89,18 @@ public class SingularSpectrumTransformTest { double x = Double.parseDouble(line); sst.update(x, outScores); printf("%f %f%n", x, outScores[0]); - if (outScores[0] > 0.95d) { + if (outScores[0] > threshold) { numChangepoints++; } } - Assert.assertTrue("#changepoints SHOULD be greater than 0: " + numChangepoints, - numChangepoints > 0); - Assert.assertTrue("#changepoints SHOULD be less than 5: " + numChangepoints, - numChangepoints < 5); + + return numChangepoints; } - @Test - public void testTwitterData() throws IOException, HiveException { + private static int detectTwitterData(@Nonnull final ScoreFunction scoreFunc, + @Nonnull final double threshold) throws IOException, HiveException { Parameters params = new Parameters(); + params.set(scoreFunc); PrimitiveObjectInspector oi = PrimitiveObjectInspectorFactory.javaDoubleObjectInspector; SingularSpectrumTransform sst = new SingularSpectrumTransform(params, oi); double[] outScores = new double[1]; @@ -76,15 +113,13 @@ public class SingularSpectrumTransformTest { double x = Double.parseDouble(line); sst.update(x, outScores); printf("%d %f %f%n", i, x, outScores[0]); - if (outScores[0] > 0.005d) { + if (outScores[0] > threshold) { numChangepoints++; } i++; } - Assert.assertTrue("#changepoints SHOULD be greater than 0: " + numChangepoints, - numChangepoints > 0); - Assert.assertTrue("#changepoints SHOULD be less than 5: " + numChangepoints, - numChangepoints < 5); + + return numChangepoints; } private static void println(String msg) { http://git-wip-us.apache.org/repos/asf/incubator-hivemall/blob/998203d5/core/src/test/java/hivemall/utils/math/MatrixUtilsTest.java ---------------------------------------------------------------------- diff --git a/core/src/test/java/hivemall/utils/math/MatrixUtilsTest.java b/core/src/test/java/hivemall/utils/math/MatrixUtilsTest.java index bc960ec..b5a5e74 100644 --- a/core/src/test/java/hivemall/utils/math/MatrixUtilsTest.java +++ b/core/src/test/java/hivemall/utils/math/MatrixUtilsTest.java @@ -19,6 +19,7 @@ package hivemall.utils.math; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.SingularValueDecomposition; import org.junit.Assert; import org.junit.Test; @@ -205,4 +206,70 @@ public class MatrixUtilsTest { Assert.assertArrayEquals(expected, actual); } + @Test + public void testPower1() { + RealMatrix A = new Array2DRowRealMatrix(new double[][] {new double[] {1, 2, 3}, new double[] {4, 5, 6}}); + + double[] x = new double[3]; + x[0] = Math.random(); + x[1] = Math.random(); + x[2] = Math.random(); + + double[] u = new double[2]; + double[] v = new double[3]; + + double s = MatrixUtils.power1(A, x, 2, u, v); + + SingularValueDecomposition svdA = new SingularValueDecomposition(A); + + Assert.assertArrayEquals(svdA.getU().getColumn(0), u, 0.001d); + Assert.assertArrayEquals(svdA.getV().getColumn(0), v, 0.001d); + Assert.assertEquals(svdA.getSingularValues()[0], s, 0.001d); + } + + @Test + public void testLanczosTridiagonalization() { + // Symmetric matrix + RealMatrix C = new Array2DRowRealMatrix(new double[][] { + new double[] {1, 2, 3, 4}, new double[] {2, 1, 4, 3}, + new double[] {3, 4, 1, 2}, new double[] {4, 3, 2, 1}}); + + // naive initial vector + double[] a = new double[] {1, 1, 1, 1}; + + RealMatrix actual = new Array2DRowRealMatrix(new double[4][4]); + MatrixUtils.lanczosTridiagonalization(C, a, actual); + + RealMatrix expected = new Array2DRowRealMatrix(new double[][] { + new double[] {40, 60, 0, 0}, new double[] {60, 10, 120, 0}, + new double[] {0, 120, 10, 120}, new double[] {0, 0, 120, 10}}); + + Assert.assertEquals(expected, actual); + } + + @Test + public void testTridiagonalEigen() { + // Tridiagonal Matrix + RealMatrix T = new Array2DRowRealMatrix(new double[][] { + new double[] {40, 60, 0, 0}, new double[] {60, 10, 120, 0}, + new double[] {0, 120, 10, 120}, new double[] {0, 0, 120, 10}}); + + double[] eigvals = new double[4]; + RealMatrix eigvecs = new Array2DRowRealMatrix(new double[4][4]); + + MatrixUtils.tridiagonalEigen(T, 2, eigvals, eigvecs); + + RealMatrix actual = eigvecs.multiply(eigvecs.transpose()); + + RealMatrix expected = new Array2DRowRealMatrix(new double[4][4]); + for (int i = 0; i < 4; i++) { + expected.setEntry(i, i, 1); + } + + for (int i = 0; i < 4; i++) { + for (int j = 0; j < 4; j++) { + Assert.assertEquals(expected.getEntry(i, j), actual.getEntry(i, j), 0.001d); + } + } + } }
