Vector previousVector = new DenseVector(currentVector.size());
Matrix basis = new SparseRowMatrix(new int[]{desiredRank, corpus.numCols()});
basis.assignRow(0, currentVector);
double alpha = 0;
double beta = 0;
DoubleMatrix2D triDiag = new DenseDoubleMatrix2D(desiredRank, desiredRank);
for (int i = 1; i < desiredRank; i++) {
startTime(TimingSection.ITERATE);
Vector nextVector = isSymmetric ? corpus.times(currentVector) : corpus.timesSquared(currentVector);
log.info("{} passes through the corpus so far...", i);
calculateScaleFactor(nextVector);
nextVector.assign(new Scale(1 / scaleFactor));
nextVector.assign(previousVector, new PlusMult(-beta));
// now orthogonalize
alpha = currentVector.dot(nextVector);
nextVector.assign(currentVector, new PlusMult(-alpha));
endTime(TimingSection.ITERATE);
startTime(TimingSection.ORTHOGANLIZE);
orthoganalizeAgainstAllButLast(nextVector, basis);
endTime(TimingSection.ORTHOGANLIZE);
// and normalize
beta = nextVector.norm(2);
if (outOfRange(beta) || outOfRange(alpha)) {
log.warn("Lanczos parameters out of range: alpha = {}, beta = {}. Bailing out early!", alpha, beta);
break;
}
final double b = beta;
nextVector.assign(new Scale(1 / b));
basis.assignRow(i, nextVector);
previousVector = currentVector;
currentVector = nextVector;
// save the projections and norms!
triDiag.set(i - 1, i - 1, alpha);
if (i < desiredRank - 1) {
triDiag.set(i - 1, i, beta);
triDiag.set(i, i - 1, beta);
}
}
startTime(TimingSection.TRIDIAG_DECOMP);
log.info("Lanczos iteration complete - now to diagonalize the tri-diagonal auxiliary matrix.");
// at this point, have tridiag all filled out, and basis is all filled out, and orthonormalized
EigenvalueDecomposition decomp = new EigenvalueDecomposition(triDiag);
DoubleMatrix2D eigenVects = decomp.getV();
DoubleMatrix1D eigenVals = decomp.getRealEigenvalues();
endTime(TimingSection.TRIDIAG_DECOMP);
startTime(TimingSection.FINAL_EIGEN_CREATE);
for (int i = 0; i < basis.numRows() - 1; i++) {
Vector realEigen = new DenseVector(corpus.numCols());
// the eigenvectors live as columns of V, in reverse order. Weird but true.
DoubleMatrix1D ejCol = eigenVects.viewColumn(basis.numRows() - i - 1);
for (int j = 0; j < ejCol.size(); j++) {
double d = ejCol.getQuick(j);
realEigen.assign(basis.getRow(j), new PlusMult(d));
}
realEigen = realEigen.normalize();