Start line:  
End line:  

Snippet Preview

Snippet HTML Code

Stack Overflow Questions
  /*
   * 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.mahout.math;
 
Cholesky decomposition shamelessly ported from JAMA.

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 final PivotedMatrix L;
   private boolean isPositiveDefinite;
 
   public CholeskyDecomposition(Matrix a) {
     this(atrue);
   }
 
   public CholeskyDecomposition(Matrix aboolean pivot) {
     int rows = a.rowSize();
      = new PivotedMatrix(new DenseMatrix(rowsrows));
 
     // must be square
     Preconditions.checkArgument(rows == a.columnSize());
 
     if (pivot) {
       decomposeWithPivoting(a);
     } else {
       decompose(a);
     }
   }
 
   private void decomposeWithPivoting(Matrix a) {
     int n = a.rowSize();
     .assign(a);
 
     // pivoted column-wise submatrix cholesky with simple pivoting
     double uberMax = .viewDiagonal().aggregate(..);
     for (int k = 0; k < nk++) {
       double max = 0;
       int pivot = k;
       for (int j = kj < nj++) {
         if (.get(jj) > max) {
           max = .get(jj);
           pivot = j;
           if (uberMax < Math.abs(max)) {
             uberMax = Math.abs(max);
           }
         }
       }
       .swap(kpivot);
 
       double akk = .get(kk);
       double epsilon = 1.0e-10 * Math.max(uberMax.viewColumn(k).aggregate(..));
 
       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
         .viewColumn(k).assign(0);
          = false;
 
         // no need to subtract from remaining sub-matrix
       } else {
         // normalize column by diagonal element
         akk = Math.sqrt(Math.max(0, akk));
         .viewColumn(k).viewPart(kn - k).assign(Functions.div(akk));
         .viewColumn(k).viewPart(0, k).assign(0);
 
         // subtract off scaled version of this column to the right
         for (int j = k + 1; j < nj++) {
           Vector columnJ = .viewColumn(j).viewPart(kn - k);
           Vector columnK = .viewColumn(k).viewPart(kn - k);
           columnJ.assign(columnK, Functions.minusMult(columnK.get(j - k)));
         }
 
       }
     }
   }
  private void decompose(Matrix a) {
    int n = a.rowSize();
    .assign(a);
    // column-wise submatrix cholesky with simple pivoting
    for (int k = 0; k < nk++) {
      double akk = .get(kk);
      // set upper part of column to 0.
      .viewColumn(k).viewPart(0, k).assign(0);
      double epsilon = 1.0e-10 * .viewColumn(k).aggregate(..);
      if (akk <= epsilon) {
        // degenerate column case.  Set diagonal to 1, all others to zero
        .viewColumn(k).viewPart(kn - k).assign(0);
         = false;
        // no need to subtract from remaining sub-matrix
      } else {
        // normalize column by diagonal element
        akk = Math.sqrt(Math.max(0, akk));
        .set(kkakk);
        .viewColumn(k).viewPart(k + 1, n - k - 1).assign(Functions.div(akk));
        // now subtract scaled version of column
        for (int j = k + 1; j < nj++) {
          Vector columnJ = .viewColumn(j).viewPart(jn - j);
          Vector columnK = .viewColumn(k).viewPart(jn - j);
          columnJ.assign(columnK, Functions.minusMult(.get(jk)));
        }
      }
    }
  }
  public boolean isPositiveDefinite() {
    return ;
  }
  public Matrix getL() {
    return .getBase();
  }
  public PivotedMatrix getPermutedL() {
    return ;
  }

  

Returns:
Returns the permutation of rows and columns that was applied to L
  public int[] getPivot() {
    return .getRowPivot();
  }
  public int[] getInversePivot() {
    return .getInverseRowPivot();
  }

  
Compute inv(L) * z efficiently.

Parameters:
z
  public Matrix solveLeft(Matrix z) {
    int n = .columnSize();
    int nx = z.columnSize();
    Matrix X = new DenseMatrix(nz.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 < ninternalK++) {
      int k = .rowUnpivot(internalK);
      for (int j = 0; j < nxj++) {
        for (int internalI = 0; internalI < internalKinternalI++) {
          int i = .rowUnpivot(internalI);
          X.set(kjX.get(kj) - X.get(ij) * .get(ki));
        }
        if (.get(kk) != 0) {
          X.set(kjX.get(kj) / .get(kk));
        } else {
          X.set(kj, 0);
        }
      }
    }
    return X;
  }

  
Compute z * inv(L') efficiently
  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 < ninternalK++) {
      int k = .rowUnpivot(internalK);
      for (int j = 0; j < nxj++) {
        for (int internalI = 0; internalI < kinternalI++) {
          int i = .rowUnpivot(internalI);
          x.set(jkx.get(jk) - x.get(ji) * .get(ki));
          if (Double.isInfinite(x.get(jk)) || Double.isNaN(x.get(jk))) {
            throw new IllegalStateException(
                String.format("Invalid value found at %d,%d (should not be possible)"jk));
          }
        }
        if (.get(kk) != 0) {
          x.set(jkx.get(jk) / .get(kk));
        } else {
          x.set(jk, 0);
        }
        if (Double.isInfinite(x.get(jk)) || Double.isNaN(x.get(jk))) {
          throw new IllegalStateException(String.format("Invalid value found at %d,%d (should not be possible)"jk));
        }
      }
    }
    return x;
  }
New to GrepCode? Check out our FAQ X