(I also sent this to the eigen@lists.tuxfamily.org list but it seems to be taking a while to appear - also there were some confusing typos in that version. Apologies for cross-posting.)

The documentation for the rankUpdate method of SparseSelfAdjointView states

Code: Select all

   /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
     * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
     *
     * \returns a reference to \c *this
     *
     * Note that it is faster to set alpha=0 than initializing the matrix to zero
     * and then keep the default value alpha=1.
     *
     * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
     * call this function with u.adjoint().
     */


The second-last paragraph is what has confused me. I want to evaluate
the Lower triangle of Z.adjoint() * Z so I am doing it as

Code: Select all

   SparseMatrix<double>  ZtZ;
   ZtZ.selfadjointView<Lower>().rankUpdate(Z.adjoint(),0.).rankUpdate(Z.adjoint());


Is that what is meant - that I should do the
.rankUpdate(Z.adjoint(), 0) then do the rankUpdate(Z.adjoint())?

A second question, if I may: What would be an effective way of
evaluating, as a sparse matrix, (I + Z'Z)? I think I would want to
initialize the result to the identity then apply the rankUpdate method
but I want to be sure that I don't accidentally create a dense, square
matrix along the way.

Note that it is faster to set alpha=0 than initializing the matrix to zero
and then keep the default value alpha=1.

hm, well... to be honest this sentence puzzle me too! It probably comes from the dense version, but even for that case, I'm not sure it makes sense. In your case simply do:

ZtZ.setZero(); // a noop if ZtZ has been just created
ZtZ.selfadjointView<Lower>().rankUpdate(Z.adjoint());

For I + Z'Z,I would rather compute Z'Z first, and then add ones to the diagonal. Unfortunately, currently for this later operation you will have to write your own loop, e.g.:

for(int j=0; j<ZtZ.outerSize(); ++j)
{
SparseMatrix<double>::InnerIterator it(ZtZ,j);
assert(it.index()==j); // because ZtZ is column major, and only the lower tri part has been filled
it.value() += 1.;
}

I'll soon add an iterator/view over the diagonal coefficients of a sparse matrix.

ggael wrote:
ZtZ.setZero(); // a noop if ZtZ has been just created
ZtZ.selfadjointView<Lower>().rankUpdate(Z.adjoint());

Interestingly, that version produces a run-time error of the form

Eigen::CwiseBinaryOp<BinaryOp, Lhs, Rhs>::CwiseBinaryOp(const Lhs&, const Rhs&, const BinaryOp&)
[with BinaryOp = Eigen::internal::scalar_sum_op<double>,
Lhs = const Eigen::SparseMatrix<double>,
Rhs = const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>,
const Eigen::SparseTriangularView<Eigen::SparseMatrix<double>, 1> >]:
Assertion `lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()' failed.

whereas the version using .rankUpdate(Z.adjoint(), 0.).rankUpdate(Z.adjoint()) succeeds. I can provide a test case if you wish.

sorry, you indeed have to resize it before since a rank update op means ZtZ = ZtZ + alpha * Z.adjoint() * Z;

that's weird the version with alpha=0 does work!!

EDIT: ok, I checked the implementation, and when alpha==0, it does:

ZtZ = Z.adjoint() * Z;

which does not make sense!! I'll fix it.

ggael wrote:

For I + Z'Z,I would rather compute Z'Z first, and then add ones to the diagonal. Unfortunately, currently for this later operation you will have to write your own loop, e.g.:

for(int j=0; j<ZtZ.outerSize(); ++j)
{
SparseMatrix<double>::InnerIterator it(ZtZ,j);
assert(it.index()==j); // because ZtZ is column major, and only the lower tri part has been filled
it.value() += 1.;
}

I'll soon add an iterator/view over the diagonal coefficients of a sparse matrix.

I think that it needs to be it.valueRef(), not it.value(), because it.value() is const Scalar&

indeed, it's valueRef()!!

dmbates wrote:

ggael wrote:

For I + Z'Z,I would rather compute Z'Z first, and then add ones to the diagonal. Unfortunately, currently for this later operation you will have to write your own loop, e.g.:

for(int j=0; j<ZtZ.outerSize(); ++j)
{
SparseMatrix<double>::InnerIterator it(ZtZ,j);
assert(it.index()==j); // because ZtZ is column major, and only the lower tri part has been filled
it.value() += 1.;
}

I'll soon add an iterator/view over the diagonal coefficients of a sparse matrix.

I think that it needs to be it.valueRef(), not it.value(), because it.value() is const Scalar&

I eventually did this in the other order, which is to initialize the result to the identity and then use selfadjointView<..>().rankUpdate(Z.adjoint())

The code you kindly provided worked except when a column of Z is all zeros, as I found out in my testing.