Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

MAHOUT-1005 - Updated to respond to Robin's comments.

  • Loading branch information...
commit ddace9fafa3cf398abb479da9b0aa1976bd767df 1 parent 0c40be5
Ted Dunning authored May 02, 2012
82  math/src/main/java/org/apache/mahout/math/solver/EigenDecomposition.java
@@ -24,6 +24,7 @@
24 24
 
25 25
 import org.apache.mahout.math.DenseMatrix;
26 26
 import org.apache.mahout.math.DenseVector;
  27
+import org.apache.mahout.math.DiagonalMatrix;
27 28
 import org.apache.mahout.math.Matrix;
28 29
 import org.apache.mahout.math.Vector;
29 30
 import org.apache.mahout.math.function.Functions;
@@ -58,10 +59,34 @@
58 59
    */
59 60
   private Matrix v;
60 61
 
  62
+  /**
  63
+   * Computes the eigen decomposition of a matrix x.  If the matrix is symmetric to
  64
+   * machine precision, a symmetric decomposition is done.  Else a complex decomposition
  65
+   * is done.  See #EigenDecomposition(Matrix, boolean) for more info.
  66
+   * @param x  The matrix to decompose.
  67
+   */
61 68
   public EigenDecomposition(Matrix x) {
62 69
     this(x, isSymmetric(x));
63 70
   }
64 71
 
  72
+  /**
  73
+   * Computes the eigen decomposition of a matrix x.  If the matrix is marked as symmetric
  74
+   * a symmetric decomposition is done.  Else a complex decomposition
  75
+   * is done.
  76
+   * <p/>
  77
+   * The decomposition has two parts, D and V.  D is a block diagonal real matrix
  78
+   * in general.  For symmetric inputs, D is entirely diagonal.  The key property
  79
+   * relating these back to the original matrix x is that
  80
+   * <p/>
  81
+   *    x V = V D
  82
+   * <p/>
  83
+   * In addition,
  84
+   * <p/>
  85
+   *    V' V = I
  86
+   * <p/>
  87
+   * @param x  The matrix to decompose.
  88
+   * @param isSymmetric true if x is to be treated as if symmetric.
  89
+   */
65 90
   public EigenDecomposition(Matrix x, boolean isSymmetric) {
66 91
     n = x.columnSize();
67 92
     d = new DenseVector(n);
@@ -94,37 +119,50 @@ public Matrix getV() {
94 119
   }
95 120
 
96 121
   /**
97  
-   * Return the real parts of the eigenvalues
  122
+   * Returns the real parts of the eigenvalues.  If the original matrix was
  123
+   * symmetric, then all eigenvalues are real and these values are identical
  124
+   * to the diagonal of the matrix returned by {@link #getD()}.
98 125
    */
99 126
   public Vector getRealEigenvalues() {
100 127
     return d;
101 128
   }
102 129
 
103 130
   /**
104  
-   * Return the imaginary parts of the eigenvalues
  131
+   * Returns the imaginary parts of the eigenvalues.  If the original matrix was
  132
+   * symmetric, then all eigenvalues are real and these values will be identically
  133
+   * zero.
105 134
    */
106 135
   public Vector getImagEigenvalues() {
107 136
     return e;
108 137
   }
109 138
 
110 139
   /**
111  
-   * Return the block diagonal eigenvalue matrix
  140
+   * Return the block diagonal eigenvalue matrix.  The eigenvalue matrix will
  141
+   * contain the real part of the eigenvalues on the diagonal and for all
  142
+   * eigenvalues that are completely real will have no off-diagonal values.
  143
+   *
  144
+   * Eigenvalues that are complex conjugate pairs will have a block of non-zero
  145
+   * values.
112 146
    *
113 147
    * @return D
114 148
    */
115 149
   public Matrix getD() {
116  
-    Matrix X = new DenseMatrix(n, n);
117  
-    X.assign(0);
118  
-    X.viewDiagonal().assign(d);
119  
-    for (int i = 0; i < n; i++) {
120  
-      final double v = e.getQuick(i);
121  
-      if (v > 0) {
122  
-        X.setQuick(i, i + 1, v);
123  
-      } else if (v < 0) {
124  
-        X.setQuick(i, i - 1, v);
  150
+    if (e.norm(1) > 0) {
  151
+      Matrix X = new DenseMatrix(n, n);
  152
+      X.assign(0);
  153
+      X.viewDiagonal().assign(d);
  154
+      for (int i = 0; i < n; i++) {
  155
+        final double v = e.getQuick(i);
  156
+        if (v > 0) {
  157
+          X.setQuick(i, i + 1, v);
  158
+        } else if (v < 0) {
  159
+          X.setQuick(i, i - 1, v);
  160
+        }
125 161
       }
  162
+      return X;
  163
+    } else {
  164
+      return new DiagonalMatrix(d);
126 165
     }
127  
-    return X;
128 166
   }
129 167
 
130 168
   // Symmetric Householder reduction to tridiagonal form.
@@ -250,7 +288,7 @@ private void tql2() {
250 288
 
251 289
     double f = 0.0;
252 290
     double tst1 = 0.0;
253  
-    double eps = Math.pow(2.0, -52.0);
  291
+    final double eps = Math.pow(2.0, -52.0);
254 292
     for (int l = 0; l < n; l++) {
255 293
 
256 294
       // Find small subdiagonal element
@@ -372,14 +410,16 @@ private Matrix orthes(Matrix x) {
372 410
 
373 411
       // Scale column.
374 412
 
375  
-      double scale = H.viewColumn(m - 1).viewPart(m, high - m + 1).norm(1);
  413
+      final Vector hColumn = H.viewColumn(m - 1).viewPart(m, high - m + 1);
  414
+      double scale = hColumn.norm(1);
376 415
 
377 416
       if (scale != 0.0) {
378 417
         // Compute Householder transformation.
379 418
 
380  
-        ort.viewPart(m, high - m + 1).assign(H.viewColumn(m - 1).viewPart(m, high - m + 1), Functions.plusMult(1 / scale));
381  
-        double h = ort.viewPart(m, high - m + 1).getLengthSquared();
  419
+        Vector ortPiece = ort.viewPart(m, high - m + 1);
  420
+        ortPiece.assign(hColumn, Functions.plusMult(1 / scale));
382 421
 
  422
+        double h = ortPiece.getLengthSquared();
383 423
         double g = Math.sqrt(h);
384 424
         if (ort.getQuick(m) > 0) {
385 425
           g = -g;
@@ -390,7 +430,6 @@ private Matrix orthes(Matrix x) {
390 430
         // Apply Householder similarity transformation
391 431
         // H = (I-u*u'/h)*H*(I-u*u')/h)
392 432
 
393  
-        Vector ortPiece = ort.viewPart(m, high - m + 1);
394 433
         for (int j = m; j < n; j++) {
395 434
           double f = ortPiece.dot(H.viewColumn(j).viewPart(m, high - m + 1)) / h;
396 435
           H.viewColumn(j).viewPart(m, high - m + 1).assign(ortPiece, Functions.plusMult(-f));
@@ -414,10 +453,11 @@ private Matrix orthes(Matrix x) {
414 453
       if (H.getQuick(m, m - 1) != 0.0) {
415 454
         ort.viewPart(m + 1, high - m).assign(H.viewColumn(m - 1).viewPart(m + 1, high - m));
416 455
         for (int j = m; j <= high; j++) {
417  
-          double g = ort.viewPart(m, high - m + 1).dot(v.viewColumn(j).viewPart(m, high - m + 1));
  456
+          final Vector ortPiece = ort.viewPart(m, high - m + 1);
  457
+          double g = ortPiece.dot(v.viewColumn(j).viewPart(m, high - m + 1));
418 458
           // Double division avoids possible underflow
419 459
           g = (g / ort.getQuick(m)) / H.getQuick(m, m - 1);
420  
-          v.viewColumn(j).viewPart(m, high - m + 1).assign(ort.viewPart(m, high - m + 1), Functions.plusMult(g));
  460
+          v.viewColumn(j).viewPart(m, high - m + 1).assign(ortPiece, Functions.plusMult(g));
421 461
         }
422 462
       }
423 463
     }
@@ -876,8 +916,6 @@ private void hqr2(Matrix h) {
876 916
     }
877 917
   }
878 918
 
879  
-
880  
-
881 919
   private static boolean isSymmetric(Matrix a) {
882 920
     /*
883 921
     Symmetry flag.

0 notes on commit ddace9f

Please sign in to comment.
Something went wrong with that request. Please try again.