How to efficiently multiply a vector by a transposed sparse matrix in Eigen?

Question

I have with a sparse matrix A and a vector x in Eigen, and I need to perform the following transposed matrix ×ばつ matrix ×ばつ vector operation: A^TAx. This can be decomposed into two matrix-vector multiplications: y=Ax followed by A^Ty. The first one is simple, but I don't know how to do the second one without the explicit construction of A^T. I can write something as:

y = A.transpose() * A * x;

However, it seems that A.transpose() explicitly constructs A^T in memory, which would not be efficient for large matrices.

Another option might be to work with a sparse A^TA, but in our case, this is a much more dense matrix than A itself.

Question: Is there any way how to multiply a vector by a transposed sparse matrix in Eigen without explicit construction of A^T or A^TA?

• Homer512: What gave you the idea that A.transpose() creates a new matrix? Most Eigen operations return an expression object or view. Specifically, it will an object of type Eigen::Transpose<Eigen::SparseMatrix<...>>. The matrix multiplication will specialize on this input parameter
• chtz: The problem is that the expression is evaluated as (A.transpose() * A) * x. Try writing A.transpose() * (A * x).
• Daniel Langr: @Homer512 IMO, Eigen is sometimes relatively poorly documented. I even couldn't find the transpose member function mentioned for the SparseMatrix class here (including the inherited member functions). However, here, it is said that: "Note that the transposition change the storage order." This is why I thought that it creates a new matrix with a different storage scheme.
• Homer512: Yes, A.transpose() * A will generate a new matrix (not immediately but by the time the matrix-vector product is evaluated). However, I understood you saying that you instead want y=Ax followed by A'y, meaning y = A.transpose() * (A * x). That only generates a temporary vector as the result of A * x, not a temporary sparse matrix. Right?
• Homer512: The generator for the documentation seems to struggle with Eigen's heavy use of templates, (partial) template specializations and the curiously recurring template pattern. At least that's my guess. My venerable old Eclipse IDE has better luck with it. But in this instance, trying to find the type, I made an educated guess, looked at the documentation for SparseMatrixBase, scrolled down to the listed header, looked at that and searched for "transpose". So yeah, the struggle is real.

Answer 1

The way I approached this was to use Visual Studio's performance profiler and have it monitor the process memory. I created a sparse square matrix of random values 25% dense. I increased the size of A from 256 to 8192 doubling it each time. The results were consistent. If A.tranpose() * (A * X) is used, the memory is not increased during that step.

Here are some results for the 8192 case

Point in Code	Memory Usage
Creating A	780 MB
A Created	1.1 GB (+333 MB)
X Created	1.1 GB (+64 KB)
Y Created	1.1 GB (+64 KB)

Since X and Y are 8192 elements I expect about a 64 KB increase for them. If A.transpose() * (A * X) created a temporary I would expect a large (400 MB ish) increase. The graph of memory is basically flat after A is created. Also the time between A to X and X to Y is roughly the same leading me to conclude there wasn't a large allocation.

The code I profiled was basically:

 Eigen::SparseMatrix<double> A(rows, cols);
 A.setFromTriplets(triplets.begin(), triplets.end());
 Eigen::VectorXd X = Eigen::Vector<double, rows>::Ones();
 Eigen::VectorXd Y = A.transpose() * (A * X);
 std::cout << Y(1) << ", " << Y(20) << '\n';