New user self-registration is currently disabled. Please email eigen-core-team @ if you need an account.
Bug 707 - In-Place decompositions
In-Place decompositions
Product: Eigen
Classification: Unclassified
Component: General
3.4 (development)
All All
: Normal Feature Request
Assigned To: Nobody
Depends on:
Blocks: 62 3.3 789
  Show dependency treegraph
Reported: 2013-11-22 13:11 UTC by Christoph Hertzberg
Modified: 2016-07-04 15:25 UTC (History)
3 users (show)

Patch for inplace decomposition via SolverName<Ref<Matrix>> (16.52 KB, patch)
2016-06-23 16:32 UTC, Gael Guennebaud
no flags Details | Diff
Usage example (1.68 KB, text/x-csrc)
2016-06-23 16:33 UTC, Gael Guennebaud
no flags Details

Description Christoph Hertzberg 2013-11-22 13:11:32 UTC
Most (dense) decompositions copy the to-be-decomposed matrix into an internal variable and decompose this variable (usually using an in-place decomposition).
Often it would be nice to expose the in-place decomposition to the public interface, to save unnecessary copies, if the original matrix is not need afterwards. For decompositions not depending on additional data (such as permutations) this should be more or less straight-forward by providing a static `computeInPlace` function.

If additional data is needed, maybe something like this could be implemented:
  Decomposition<Ref<MatrixXd> > decomp(mat);
Comment 1 Gael Guennebaud 2014-09-08 11:02:50 UTC
The first issue to fix is whether we want to share the *input* matrix only (or matrices, e.g., for general eigen problems) or all *outputs*?

Since the main use case is to save memory and copies when the input is not needed anymore, it seems to me that sharing inputs only would be enough.

Regarding the implementation, we could add a computeInPlace() non-static member function. Internally, the decomposition object would possesses both a plain Matrix<> m_mat object and a Ref<> m_ref object. If compute() is called, then the m_mat will be resized and m_ref will reference it. If computeInPlace is called, the m_ref will reference the input matrix.

The main advantage of this strategy is its flexibility as the same decomposition object can be used for both usages. On the other hand, this strategy does not permit to save memory for fixed size objects (no big deal), and the options of the  Ref<> object have to be fixed by us. They include:
- Alignment (default OFF)
- OuterStride (default Dynamic)
- InnerStride (default 1)

With your suggestions, i.e., Decomposition<Ref<MatrixXd> >, then the power-user would be able to customize these options but the "inplace" nature of a call to compute() would be less obvious.
Comment 2 Christoph Hertzberg 2015-06-16 19:56:25 UTC
This will actually be necessary if we want an allocation-free LAPACK implementation (we'd also need to pass permutation vectors by reference for that).
Comment 3 Gael Guennebaud 2015-06-16 21:38:45 UTC
Actually, this is already the case for Cholesky and LU for which our LAPACK interface uses our internal in-place routines and completely ignore the respective public classes.
Comment 4 Rasmus Munk Larsen 2016-03-22 18:09:53 UTC
Bump. This is becoming a major pain point for certain memory constrained problems that I work on.
Comment 5 Gael Guennebaud 2016-06-14 15:19:10 UTC
As discussed there:,
declaring PartialPivLU<Ref<MatrixXd> > almost already works. We only have to add a non const ctor like:

template<typename MatrixType>
template<typename InputType>
PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix)
  : m_lu(matrix.derived()),

Of course, this does not cover the other objects like the permutations and workspace, but those are usually much less concerning.

Another limitation is that you cannot call compute() with another matrix to work inplace in the new one:

MatrixXd A, B;

PartialPivLU<Ref<MatrixXd> > lu(A); // works inplace in A
lu.compute(B); // copy B in A and works in A

Actually, it is not clear what this second line should do, and copying in A could be the desired behavior. We could address this issue by adding a computeInPlace() method, which would be enabled only if the MatrixType is a Ref<>.
Comment 6 Gael Guennebaud 2016-06-23 16:32:39 UTC
Created attachment 714 [details]
Patch for inplace decomposition via SolverName<Ref<Matrix>>

This patch adds support for inplace decomposition for: all LU, Cholesky, and QR decompositions.
Comment 7 Gael Guennebaud 2016-06-23 16:33:40 UTC
Created attachment 715 [details]
Usage example

Here is a simple exemple checking and demonstrating inplace decomposition.
Comment 8 Gael Guennebaud 2016-06-23 16:34:35 UTC
Since the required changes are so small, I'd like to add it for 3.3.
Comment 9 Christoph Hertzberg 2016-06-23 17:07:24 UTC
>     template<typename InputType>
>     explicit LDLT(const EigenBase<InputType>& matrix)
>       : m_matrix(matrix.rows(), matrix.cols()),

Can't we just replace the line above by m_matrix(matrix.derived()), and save the additional constructor?

As this would cause a second copy (unless the compiler can optimize it away), how about adding another (protected/private) compute() method without parameters, and let the current compute(const EigenBase<...>&) method just make a copy and call that new method?

This would also allow calling in-place decompositions on blocks of matrices or temporarily created Maps.

The same goes for the other decompositions, of course.
Comment 10 Gael Guennebaud 2016-06-23 20:13:42 UTC
The "problem" is that we need a non-const ctor for the inplace version.

The extra copy already exists because the non-const ctor can still be called even without a Ref<> as MatrixType, so yes we should fix it and your suggestion is a good one.
Comment 11 Gael Guennebaud 2016-07-04 15:25:21 UTC
Summary:     Bug 707: add inplace decomposition through Ref<> for Cholesky, LU and QR decompositions.


Good enough for me.

If someone needs more control, such as sharing storage of permutations and Householder coefficients, then let's open a new bug entry.

Note You need to log in before you can comment on or make changes to this bug.