You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@hivemall.apache.org by my...@apache.org on 2016/12/02 08:02:20 UTC

[30/50] [abbrv] incubator-hivemall git commit: Support implicit-Krylov-approximation-based efficient SST

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 <k....@gmail.com>
Authored: Wed Sep 28 19:49:48 2016 +0900
Committer: Takuya Kitazawa <k....@gmail.com>
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);
+            }
+        }
+    }
 }