in src/main/java/org/apache/datasketches/vector/decomposition/MatrixOpsImplOjAlgo.java [174:209]
private void computeSISVD(final MatrixStore<Double> A, final boolean computeVectors) {
// want to iterate on smaller dimension of A (n x d)
// currently, error in constructor if d < n, so n is always the smaller dimension
if (block_ == null) {
block_ = Primitive64Store.FACTORY.makeFilled(d_, k_, new Normal(0.0, 1.0));
qr_ = QR.PRIMITIVE.make(block_);
T_ = Primitive64Store.FACTORY.make(n_, k_);
} else {
block_.fillAll(new Normal(0.0, 1.0));
}
// orthogonalize for numeric stability
qr_.decompose(block_);
qr_.getQ().supplyTo(block_);
for (int i = 0; i < numSISVDIter_; ++i) {
A.multiply(block_, T_);
// again, just for stability
qr_.decompose(T_.premultiply(A.transpose()));
qr_.getQ().supplyTo(block_);
}
// Rayleigh-Ritz postprocessing
final SingularValue<Double> svd = SingularValue.PRIMITIVE.make(T_);
svd.compute(block_.premultiply(A));
svd.getSingularValues(sv_);
if (computeVectors) {
// V = block * Q2^T so V^T = Q2 * block^T
// and ojAlgo figures out that it only needs to fill the first k_ rows of Vt_
svd.getV().multiply(block_.transpose()).supplyTo(Vt_);
}
}