public CholeskyDecomposition regularizedCholesky(double[][] gram, int max_attempts, boolean throw_exception) { int attempts = 0; double addedL2 = 0; Matrix gmat = new Matrix(gram); CholeskyDecomposition chol = new CholeskyDecomposition(gmat); while (!chol.isSPD() && attempts < max_attempts) { if (addedL2 == 0) addedL2 = 1e-5; else addedL2 *= 10; ++attempts; for (int i = 0; i < gram.length; i++) gmat.set(i,i,addedL2); // try to add L2 penalty to make the Gram SPD Log.info("Added L2 regularization = " + addedL2 + " to diagonal of Gram matrix"); chol = new CholeskyDecomposition(gmat); } if (!chol.isSPD() && throw_exception) throw new Gram.NonSPDMatrixException(); return chol; }
Matrix tmp = _chol.solve(new Matrix(new double[][] {xrow}).transpose()); xrow = tmp.getColumnPackedCopy();
public Matrix solveSPD(Matrix b) { if (b instanceof JamaDenseDoubleMatrix2D) { JamaDenseDoubleMatrix2D b2 = (JamaDenseDoubleMatrix2D) b; CholeskyDecomposition chol = new CholeskyDecomposition(matrix); return new JamaDenseDoubleMatrix2D(chol.solve(b2.matrix)); } else { return super.solve(b); } }
public Matrix chol() { CholeskyDecomposition chol = new CholeskyDecomposition(matrix); Matrix r = new JamaDenseDoubleMatrix2D(chol.getL()); return r; }
protected void cacheValues() { inv_covar = covar.inverse(); pdf_const_factor = 1.0 / Math.sqrt((Math.pow((2 * Math.PI), N) * covar.det())); chol = covar.chol().getL(); }
private void initialXClosedForm(DataInfo dinfo, Archetypes yt_arch, double[] normSub, double[] normMul) { Log.info("Initializing X = AY'(YY' + gamma I)^(-1) where A = training data"); double[][] ygram = ArrayUtils.formGram(yt_arch._archetypes); if (_parms._gamma_y > 0) { for (int i = 0; i < ygram.length; i++) ygram[i][i] += _parms._gamma_y; } CholeskyDecomposition yychol = regularizedCholesky(ygram, 10, false); if(!yychol.isSPD()) Log.warn("Initialization failed: (YY' + gamma I) is non-SPD. Setting initial X to standard normal" + " random matrix. Results will be numerically unstable"); else { CholMulTask cmtsk = new CholMulTask(yychol, yt_arch, _ncolA, _ncolX, dinfo._cats, normSub, normMul); cmtsk.doAll(dinfo._adaptedFrame); } }
/** Cholesky Decomposition @return CholeskyDecomposition @see CholeskyDecomposition */ public CholeskyDecomposition chol () { return new CholeskyDecomposition(this); }
public Matrix solveSPD(Matrix b) { if (b instanceof JamaDenseDoubleMatrix2D) { JamaDenseDoubleMatrix2D b2 = (JamaDenseDoubleMatrix2D) b; CholeskyDecomposition chol = new CholeskyDecomposition(matrix); return new JamaDenseDoubleMatrix2D(chol.solve(b2.matrix)); } else { return super.solve(b); } }
public Matrix chol() { CholeskyDecomposition chol = new CholeskyDecomposition(matrix); Matrix r = new JamaDenseDoubleMatrix2D(chol.getL()); return r; }
@Override public double[] sample(Random rng) { final int N = mean.getColumnDimension(); final Matrix chol = getCovariance().chol().getL(); final Matrix vec = new Matrix(N, 1); for (int i = 0; i < N; i++) vec.set(i, 0, rng.nextGaussian()); final Matrix result = this.mean.plus(chol.times(vec).transpose()); return result.getArray()[0]; }
private void initialXClosedForm(DataInfo dinfo, Archetypes yt_arch, double[] normSub, double[] normMul) { Log.info("Initializing X = AY'(YY' + gamma I)^(-1) where A = training data"); double[][] ygram = ArrayUtils.formGram(yt_arch._archetypes); if (_parms._gamma_y > 0) { for (int i = 0; i < ygram.length; i++) ygram[i][i] += _parms._gamma_y; } CholeskyDecomposition yychol = regularizedCholesky(ygram, 10, false); if(!yychol.isSPD()) Log.warn("Initialization failed: (YY' + gamma I) is non-SPD. Setting initial X to standard normal" + " random matrix. Results will be numerically unstable"); else { CholMulTask cmtsk = new CholMulTask(yychol, yt_arch, _ncolA, _ncolX, dinfo._cats, normSub, normMul); cmtsk.doAll(dinfo._adaptedFrame); } }
public static CholeskyDecomposition cholesky(double[][] v) { return new CholeskyDecomposition(new Matrix(v)); }
public void test () { Log.info("CholTest::test enter"); for (int sz = 6000; sz < 10000; sz+=2000) { Log.info("CholTest::test sz is " + sz); DataSetup data = new DataSetup(sz, 12345); long start = System.currentTimeMillis(); CholeskyDecomposition jamaChol = new Matrix(data.xx).chol(); Log.info("JAMA CHOLESKY [N = " + sz + "] TAKES " + (System.currentTimeMillis() - start) + " MILLISECONDS."); if (!jamaChol.isSPD()) continue; ForkJoinPool fjp = new ForkJoinPool(32); for (int t = 2; t <= 32; t += 2) { for (int step : STEPS) fjp.invoke(new TestSetup(new DataSetup(data.xx),jamaChol.getL().getArray(),step,t)); } } Log.info("CholTest::test exit"); }
public Matrix invSPD() { CholeskyDecomposition chol = new CholeskyDecomposition(matrix); return new JamaDenseDoubleMatrix2D(chol.solve(Jama.Matrix.identity(matrix.getRowDimension(), matrix.getRowDimension()))); }
public CholeskyDecomposition regularizedCholesky(double[][] gram, int max_attempts, boolean throw_exception) { int attempts = 0; double addedL2 = 0; Matrix gmat = new Matrix(gram); CholeskyDecomposition chol = new CholeskyDecomposition(gmat); while (!chol.isSPD() && attempts < max_attempts) { if (addedL2 == 0) addedL2 = 1e-5; else addedL2 *= 10; ++attempts; for (int i = 0; i < gram.length; i++) gmat.set(i,i,addedL2); // try to add L2 penalty to make the Gram SPD Log.info("Added L2 regularization = " + addedL2 + " to diagonal of Gram matrix"); chol = new CholeskyDecomposition(gmat); } if (!chol.isSPD() && throw_exception) throw new Gram.NonSPDMatrixException(); return chol; }
public double [] sample(Random random) { Matrix L = getChol().getL(); // start with a vector of iid std normals Matrix stdNormal = new Matrix(dim(), 1); for (int i = 0; i < dim(); i++) { stdNormal.set(i, 0, SampleUtils.sampleGaussian(random)); } Matrix result = L.times(stdNormal); result.plusEquals(mean); return result.getColumnPackedCopy(); } public double[] sampleObject(Random random) { return sample(random); }
@Override public long process(BenchmarkMatrix[] inputs, BenchmarkMatrix[] outputs, long numTrials) { Matrix matA = inputs[0].getOriginal(); Matrix result = null; int N = matA.getColumnDimension(); long prev = System.nanoTime(); for( long i = 0; i < numTrials; i++ ) { result = matA.chol().solve(Matrix.identity(N,N)); } long elapsed = System.nanoTime()-prev; if( outputs != null ) { outputs[0] = new JamaBenchmarkMatrix(result); } return elapsed; } }
@Override public long process(BenchmarkMatrix[] inputs, BenchmarkMatrix[] outputs, long numTrials) { Matrix matA = inputs[0].getOriginal(); Matrix L = null; long prev = System.nanoTime(); for( long i = 0; i < numTrials; i++ ) { CholeskyDecomposition chol = matA.chol(); if( !chol.isSPD() ) { throw new DetectedException("Is not SPD"); } L = chol.getL(); } long elapsed = System.nanoTime()-prev; if( outputs != null ) { outputs[0] = new JamaBenchmarkMatrix(L); } return elapsed; } }