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.");