Author: tdunning
Date: Thu Dec 22 23:51:27 2011
New Revision: 1222513
URL: http://svn.apache.org/viewvc?rev=1222513&view=rev
Log:
MAHOUT-792 - Cholesky decomposition
MAHOUT-792 - New Cholesky decomposition for SSVD update.
MAHOUT-792 - Switch Cholesky to views.
Added:
mahout/trunk/math/src/main/java/org/apache/mahout/math/CholeskyDecomposition.java
mahout/trunk/math/src/test/java/org/apache/mahout/math/CholeskyDecompositionTest.java
Added:
mahout/trunk/math/src/main/java/org/apache/mahout/math/CholeskyDecomposition.java
URL:
http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/CholeskyDecomposition.java?rev=1222513&view=auto
==============================================================================
---
mahout/trunk/math/src/main/java/org/apache/mahout/math/CholeskyDecomposition.java
(added)
+++
mahout/trunk/math/src/main/java/org/apache/mahout/math/CholeskyDecomposition.java
Thu Dec 22 23:51:27 2011
@@ -0,0 +1,211 @@
+package org.apache.mahout.math;
+
+import com.google.common.base.Preconditions;
+import org.apache.mahout.math.function.Functions;
+
+/**
+ * Cholesky decomposition shamelessly ported from JAMA.
+ * <p/>
+ * A Cholesky decomposition of a semi-positive definite matrix A is a lower
triangular matrix L such
+ * that L L^* = A. If A is full rank, L is unique. If A is real, then it
must be symmetric and R
+ * will also be real.
+ */
+public class CholeskyDecomposition {
+ private PivotedMatrix L;
+ private boolean isPositiveDefinite;
+
+ public CholeskyDecomposition(Matrix a) {
+ this(a, true);
+ }
+
+ public CholeskyDecomposition(Matrix a, boolean pivot) {
+ int rows = a.rowSize();
+ L = new PivotedMatrix(new DenseMatrix(rows, rows));
+
+ // must be square
+ Preconditions.checkArgument(rows == a.columnSize());
+
+ if (pivot) {
+ decomposeWithPivoting(a);
+ } else {
+ decompose(a);
+ }
+ }
+
+ private void decomposeWithPivoting(Matrix a) {
+ int n = a.rowSize();
+ L.assign(a);
+
+ // pivoted column-wise submatrix cholesky with simple pivoting
+ double uberMax = L.viewDiagonal().aggregate(Functions.MAX, Functions.ABS);
+ for (int k = 0; k < n; k++) {
+ double max = 0;
+ int pivot = k;
+ for (int j = k; j < n; j++) {
+ if (L.get(j, j) > max) {
+ max = L.get(j, j);
+ pivot = j;
+ if (uberMax < Math.abs(max)) {
+ uberMax = Math.abs(max);
+ }
+ }
+ }
+ L.swap(k, pivot);
+
+ double akk = L.get(k, k);
+ double epsilon = 1e-10 * Math.max(uberMax,
L.viewColumn(k).aggregate(Functions.MAX, Functions.ABS));
+
+ if (akk < -epsilon) {
+ // can't have decidedly negative element on diagonal
+ throw new IllegalArgumentException("Matrix is not positive
semi-definite");
+ } else if (akk <= epsilon) {
+ // degenerate column case. Set all to zero
+ L.viewColumn(k).assign(0);
+ isPositiveDefinite = false;
+
+ // no need to subtract from remaining sub-matrix
+ } else {
+ // normalize column by diagonal element
+ akk = Math.sqrt(Math.max(0, akk));
+ L.viewColumn(k).viewPart(k, n - k).assign(Functions.div(akk));
+ L.viewColumn(k).viewPart(0, k).assign(0);
+
+ // subtract off scaled version of this column to the right
+ for (int j = k + 1; j < n; j++) {
+ Vector columnJ = L.viewColumn(j).viewPart(k, n - k);
+ Vector columnK = L.viewColumn(k).viewPart(k, n - k);
+ columnJ.assign(columnK, Functions.minusMult(columnK.get(j - k)));
+ }
+
+ }
+ }
+ }
+
+ private void decompose(Matrix a) {
+ int n = a.rowSize();
+ L.assign(a);
+
+ // column-wise submatrix cholesky with simple pivoting
+ for (int k = 0; k < n; k++) {
+
+ double akk = L.get(k, k);
+
+ // set upper part of column to 0.
+ L.viewColumn(k).viewPart(0, k).assign(0);
+
+ double epsilon = 1e-10 * L.viewColumn(k).aggregate(Functions.MAX,
Functions.ABS);
+ if (akk <= epsilon) {
+ // degenerate column case. Set diagonal to 1, all others to zero
+ L.viewColumn(k).viewPart(k, n - k).assign(0);
+
+ isPositiveDefinite = false;
+
+ // no need to subtract from remaining sub-matrix
+ } else {
+ // normalize column by diagonal element
+ akk = Math.sqrt(Math.max(0, akk));
+ L.set(k, k, akk);
+ L.viewColumn(k).viewPart(k + 1, n - k - 1).assign(Functions.div(akk));
+
+ // now subtract scaled version of column
+ for (int j = k + 1; j < n; j++) {
+ Vector columnJ = L.viewColumn(j).viewPart(j, n - j);
+ Vector columnK = L.viewColumn(k).viewPart(j, n - j);
+ columnJ.assign(columnK, Functions.minusMult(L.get(j, k)));
+ }
+ }
+ }
+ }
+
+ public boolean isPositiveDefinite() {
+ return isPositiveDefinite;
+ }
+
+ public Matrix getL() {
+ return L.getBase();
+ }
+
+ public PivotedMatrix getPermutedL() {
+ return L;
+ }
+
+ /**
+ * @return Returns the permutation of rows and columns that was applied to L
+ */
+ public int[] getPivot() {
+ return L.getRowPivot();
+ }
+
+ public int[] getInversePivot() {
+ return L.getInverseRowPivot();
+ }
+
+ /**
+ * Compute inv(L) * z efficiently.
+ *
+ * @param z
+ */
+ public Matrix solveLeft(Matrix z) {
+ int n = L.columnSize();
+ int nx = z.columnSize();
+
+ Matrix X = new DenseMatrix(n, z.columnSize());
+ X.assign(z);
+
+ // Solve L*Y = Z using back-substitution
+ // note that k and i have to go in a funny order because L is pivoted
+ for (int internalK = 0; internalK < n; internalK++) {
+ int k = L.rowUnpivot(internalK);
+ for (int j = 0; j < nx; j++) {
+ for (int internalI = 0; internalI < internalK; internalI++) {
+ int i = L.rowUnpivot(internalI);
+ X.set(k, j, X.get(k, j) - X.get(i, j) * L.get(k, i));
+ }
+ if (L.get(k, k) != 0) {
+ X.set(k, j, X.get(k, j) / L.get(k, k));
+ } else {
+ X.set(k, j, 0);
+ }
+ }
+ }
+ return X;
+ }
+
+ /**
+ * Compute z * inv(L') efficiently
+ *
+ * @param z
+ * @return
+ */
+ public Matrix solveRight(Matrix z) {
+ int n = z.columnSize();
+ int nx = z.rowSize();
+
+ Matrix x = new DenseMatrix(z.rowSize(), z.columnSize());
+ x.assign(z);
+
+ // Solve Y*L' = Z using back-substitution
+ for (int internalK = 0; internalK < n; internalK++) {
+ int k = L.rowUnpivot(internalK);
+ for (int j = 0; j < nx; j++) {
+ for (int internalI = 0; internalI < k; internalI++) {
+ int i = L.rowUnpivot(internalI);
+ x.set(j, k, x.get(j, k) - x.get(j, i) * L.get(k, i));
+ if (Double.isInfinite(x.get(j, k)) || Double.isNaN(x.get(j, k))) {
+ System.out.printf("bad at %d,%d\n", j, k);
+ }
+ }
+ if (L.get(k, k) != 0) {
+ x.set(j, k, x.get(j, k) / L.get(k, k));
+ } else {
+ x.set(j, k, 0);
+ }
+ if (Double.isInfinite(x.get(j, k)) || Double.isNaN(x.get(j, k))) {
+ System.out.printf("bad at %d,%d\n", j, k);
+ }
+ }
+ }
+ return x;
+ }
+}
+
Added:
mahout/trunk/math/src/test/java/org/apache/mahout/math/CholeskyDecompositionTest.java
URL:
http://svn.apache.org/viewvc/mahout/trunk/math/src/test/java/org/apache/mahout/math/CholeskyDecompositionTest.java?rev=1222513&view=auto
==============================================================================
---
mahout/trunk/math/src/test/java/org/apache/mahout/math/CholeskyDecompositionTest.java
(added)
+++
mahout/trunk/math/src/test/java/org/apache/mahout/math/CholeskyDecompositionTest.java
Thu Dec 22 23:51:27 2011
@@ -0,0 +1,139 @@
+package org.apache.mahout.math;
+
+import org.apache.mahout.common.RandomUtils;
+import org.apache.mahout.math.function.DoubleDoubleFunction;
+import org.apache.mahout.math.function.DoubleFunction;
+import org.apache.mahout.math.function.Functions;
+import org.junit.Assert;
+import org.junit.Test;
+
+import java.util.Random;
+
+public class CholeskyDecompositionTest extends MahoutTestCase {
+ @Test
+ public void rank1() {
+ Matrix x = new DenseMatrix(3, 3);
+ x.viewRow(0).assign(new double[]{1, 2, 3});
+ x.viewRow(1).assign(new double[]{2, 4, 6});
+ x.viewRow(2).assign(new double[]{3, 6, 9});
+
+ CholeskyDecomposition rr = new
CholeskyDecomposition(x.transpose().times(x), false);
+ assertEquals(0, new DenseVector(new double[]{3.741657, 7.483315,
11.22497}).aggregate(rr.getL().transpose().viewRow(0), Functions.PLUS, new
DoubleDoubleFunction() {
+ @Override
+ public double apply(double arg1, double arg2) {
+ return Math.abs(arg1) - Math.abs(arg2);
+ }
+ }), 1e-5);
+
+ assertEquals(0, rr.getL().viewPart(0, 3, 1, 2).aggregate(Functions.PLUS,
Functions.ABS), 1e-9);
+ }
+
+ @Test
+ public void test1() {
+ Matrix lastL;
+
+ final Random rand = RandomUtils.getRandom();
+
+ Matrix z = new DenseMatrix(100, 100);
+ z.assign(new DoubleFunction() {
+ @Override
+ public double apply(double arg1) {
+ return rand.nextDouble();
+ }
+ });
+
+ Matrix A = z.times(z.transpose());
+
+ for (boolean type = false; !type; type=true) {
+ CholeskyDecomposition cd = new CholeskyDecomposition(A, type);
+ Matrix L = cd.getL();
+ lastL = L;
+// Assert.assertTrue("Positive definite", cd.isPositiveDefinite());
+
+ Matrix Abar = L.times(L.transpose());
+
+ double error = A.minus(Abar).aggregate(Functions.MAX, Functions.ABS);
+ Assert.assertEquals("type = " + type, 0, error, 1e-10);
+
+ // L should give us a quick and dirty LQ decomposition
+ Matrix q = cd.solveLeft(z);
+ Matrix id = q.times(q.transpose());
+ for (int i = 0; i < id.columnSize(); i++) {
+ Assert.assertEquals("type = " + type, 1, id.get(i, i), 1e-9);
+ Assert.assertEquals("type = " + type, 1, id.viewRow(i).norm(1), 1e-9);
+ }
+
+ // and QR as well
+ q = cd.solveRight(z.transpose());
+ id = q.transpose().times(q);
+ for (int i = 0; i < id.columnSize(); i++) {
+ Assert.assertEquals("type = " + type, 1, id.get(i, i), 1e-9);
+ Assert.assertEquals("type = " + type, 1, id.viewRow(i).norm(1), 1e-9);
+ }
+ }
+ }
+
+ @Test
+ public void test2() {
+ // Test matrix from Nicholas Higham's paper at
http://eprints.ma.man.ac.uk/1199/01/covered/MIMS_ep2008_116.pdf
+ double[][] values = new double[3][];
+ values[0] = new double[]{1, -1, 1};
+ values[1] = new double[]{-1, 1, -1};
+ values[2] = new double[]{1, -1, 2};
+
+ Matrix A = new DenseMatrix(values);
+
+ // without pivoting
+ CholeskyDecomposition cd = new CholeskyDecomposition(A, false);
+ assertEquals(0,
cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS,
Functions.ABS), 1e-10);
+
+ // with pivoting
+ cd = new CholeskyDecomposition(A);
+ assertEquals(0,
cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS,
Functions.ABS), 1e-10);
+ }
+
+
+ @Test
+ public void testRankDeficient() {
+ Matrix A = rank4Matrix();
+
+ CholeskyDecomposition cd = new CholeskyDecomposition(A, false);
+
+ cd = new CholeskyDecomposition(A);
+
+ PivotedMatrix Ax = new PivotedMatrix(A, cd.getPivot());
+ CholeskyDecomposition cd2 = new CholeskyDecomposition(Ax, false);
+
+ assertEquals(0,
cd2.getL().times(cd2.getL().transpose()).minus(Ax).aggregate(Functions.PLUS,
Functions.ABS), 1e-10);
+ assertEquals(0,
cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS,
Functions.ABS), 1e-10);
+
+ Assert.assertFalse(cd.isPositiveDefinite());
+ Matrix L = cd.getL();
+ Matrix Abar = L.times(L.transpose());
+ double error = A.minus(Abar).aggregate(Functions.MAX, Functions.ABS);
+ Assert.assertEquals(0, error, 1e-10);
+ }
+
+ private Matrix rank4Matrix() {
+ final Random rand = RandomUtils.getRandom();
+
+ Matrix u = new DenseMatrix(10, 4);
+ u.assign(new DoubleFunction() {
+ @Override
+ public double apply(double arg1) {
+ return rand.nextDouble();
+ }
+ });
+
+ Matrix v = new DenseMatrix(10, 4);
+ v.assign(new DoubleFunction() {
+ @Override
+ public double apply(double arg1) {
+ return rand.nextDouble();
+ }
+ });
+
+ Matrix z = u.times(v.transpose());
+ return z.times(z.transpose());
+ }
+}
