private void computeSISVD()

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_);
    }
  }