In some applications I and my colleagues have been working on, it would be useful to be able to make low-rank updates to an existing factorization, in cases where we don't have access to an LLT or LDLT instance that was used to compute it. Currently, the only route to the rankUpdate methods are through the public interface of an LLT/LDLT instance, and the only way to construct an instance of LLT/LDLT is to instantiate it with an SP[S]D matrix, which is then factorized in the various constructors. I can see at least two routes to enabling what we want 1. Add an extra constructor to the LLT/LDLT classes, which takes a tril matrix and skips the computation 2. Refactor things to separate rankUpdate functions from the LLT/LDLT classes entirely Other options that come to mind amount to minor variations on these (e.g., instead of a new constructor, use a factory method, an additional boolean flag on existing constructors, or a separate boolean trait). I implemented a simplistic version of the "new constructor" approach in my local fork; and it is straightforward to implement and use. My current sketch of an implementation is pretty naïve, though: eg, it ignores existence of _UpLo trait, and in order to avoid collision with existing constructor signatures, introduces a boolean flag 'validate_input' that lets me tell the ctor whether to validate that the input is indeed triangular (it's not clear to me, as a total Eigen dev n00b whether such optional validations are customary in this code base; anyway, it feels like a sloppy hack to differentiate this ctor from the others). I'd love to get feedback on the above sketch of a proposal. If it sounds generally reasonable, I can actually write a less sketch-y proposal and eventually start submitting patches. I'd be happy to contribute the work here, as this will be very useful to us and, I hope, to others. Thanks in advance!
For LLT you can call: Eigen::internal::llt_rank_update_lower(L,v,s); which updates L such that: L*L' = L*L' + s * v * v' For LDLT: Eigen::internal::ldlt_inplace<Lower>::updateInPlace(L,v,s);
Ah I see, thanks Gael. Our convention in Google cpp code is that 'internal' namespaces shouldn't be used (indeed, usually not even visible) outside the module in which they're defined. Is this supported in Eigen? Should we extract these methods outside the internal namespace? My worry would be that such internal interfaces would not be protected from API change in the future.
yes, of course playing with Eigen::internal is safe regarding API stability. If we don't mind a copy, we could simply offer: MatrixXd L; LLT<MatrixXd> llt_wrapper; llt_wrapper.initFromFactors(L); For LDLT, it's a bit more tricky because there is the diagonal vector (internally stored in the diagonal!) and the permutation. We could be flexible and allows for any combination: LDLT<MatrixXd> ldlt_wrapper; ldlt_wrapper.initFromFactors(LD); // P=identity ldlt_wrapper.initFromFactors(L,D); // P=identity ldlt_wrapper.initFromFactors(LD,P); ldlt_wrapper.initFromFactors(L,D,P); We could also imagine using the "in-place" approach (https://eigen.tuxfamily.org/dox/group__InplaceDecomposition.html) to avoid copies, but in this case we need a tagged ctor (or a factory, but I4d rather lean for the ctor): MatrixXd L; LLT<Ref<MatrixXd>> llt_wrapper(InitFromFactors, L); LDLT<Ref<MatrixXd>> ldlt_wrapper(InitFromFactors, LD);
(In reply to Gael Guennebaud from comment #3) > yes, of course playing with Eigen::internal is safe regarding API stability. I guess you forgot a "not" :) > MatrixXd L; > LLT<Ref<MatrixXd>> llt_wrapper(InitFromFactors, L); > LDLT<Ref<MatrixXd>> ldlt_wrapper(InitFromFactors, LD); If we offer these constructors, we could offer the same constructors for the standard "copy" variant. But indeed, we should additionally offer the initFromFactors() methods (which would just be called by the constructors). The alternative with overloaded constructors for TriangularView<...>, LDLT<MatrixXd> ldlt_wrapper(LD.triangularView<UnitUpper>(), LD.diagonal(), P); might be a bit dangerous for generic code. A very hacky solution, would be to const_cast the matrixLLT() return value and assign (or swap) an already decomposed matrix to it: MatrixXd L(2,2); L << 1, 0, 2, 3; Vector2d x(1,2), b = L * L.transpose() * x; LLT<MatrixXd> llt( (MatrixXd()) ); // extra parentheses because C++ ... const_cast<MatrixXd&>(llt.matrixLLT())=L; std::cout << llt.solve(b) << "\n"; With `LLT<Ref<MatrixXd> >`, the const_cast could be avoided ... For the LDLT case (and other decompositions), this would require keeping the permutations vector at the same size, or strange things will happen. And needing const_cast generally is not a sign of a nice API (to say it nicely).
(In reply to Christoph Hertzberg from comment #4) > (In reply to Gael Guennebaud from comment #3) > > yes, of course playing with Eigen::internal is safe regarding API stability. > > I guess you forgot a "not" :) Thanks for clarifying -- I was not so sure how to read this! > > > > MatrixXd L; > > LLT<Ref<MatrixXd>> llt_wrapper(InitFromFactors, L); > > LDLT<Ref<MatrixXd>> ldlt_wrapper(InitFromFactors, LD); > > If we offer these constructors, we could offer the same constructors for the > standard "copy" variant. But indeed, we should additionally offer the > initFromFactors() methods (which would just be called by the constructors). I like the idea of tagged constructors (hadn't read about tag dispatch before -- nice trick :)) > > > The alternative with overloaded constructors for TriangularView<...>, > > LDLT<MatrixXd> ldlt_wrapper(LD.triangularView<UnitUpper>(), LD.diagonal(), > P); > > might be a bit dangerous for generic code. > > A very hacky solution, would be to const_cast the matrixLLT() return value > and assign (or swap) an already decomposed matrix to it: > > MatrixXd L(2,2); L << 1, 0, 2, 3; > Vector2d x(1,2), b = L * L.transpose() * x; > > LLT<MatrixXd> llt( (MatrixXd()) ); // extra parentheses because C++ ... > const_cast<MatrixXd&>(llt.matrixLLT())=L; > > std::cout << llt.solve(b) << "\n"; > > With `LLT<Ref<MatrixXd> >`, the const_cast could be avoided ... > > For the LDLT case (and other decompositions), this would require keeping the > permutations vector at the same size, or strange things will happen. > And needing const_cast generally is not a sign of a nice API (to say it > nicely). Agreed these solutions sound a bit suboptimal. One advantage of leaning on the type system in the TriangularView approach is it obviates any need for validation of the inputs. We could also avoid this by simply declaring that the implementation only looks at the lower/upper triangular portion of the input, and skips validation altogether. This seems more in the spirit of the current implementation, which only performs validation in the sense of setting m_info when a computation goes awry.
(In reply to Christopher Suter from comment #5) > (In reply to Christoph Hertzberg from comment #4) > > (In reply to Gael Guennebaud from comment #3) > > > yes, of course playing with Eigen::internal is safe regarding API stability. > > > > I guess you forgot a "not" :) > > Thanks for clarifying -- I was not so sure how to read this! > > > > > > > > MatrixXd L; > > > LLT<Ref<MatrixXd>> llt_wrapper(InitFromFactors, L); > > > LDLT<Ref<MatrixXd>> ldlt_wrapper(InitFromFactors, LD); > > > > If we offer these constructors, we could offer the same constructors for the > > standard "copy" variant. But indeed, we should additionally offer the > > initFromFactors() methods (which would just be called by the constructors). > > I like the idea of tagged constructors (hadn't read about tag dispatch > before -- nice trick :)) > > > > > > > > The alternative with overloaded constructors for TriangularView<...>, > > > > LDLT<MatrixXd> ldlt_wrapper(LD.triangularView<UnitUpper>(), LD.diagonal(), > > P); > > > > might be a bit dangerous for generic code. > > > > A very hacky solution, would be to const_cast the matrixLLT() return value > > and assign (or swap) an already decomposed matrix to it: > > > > MatrixXd L(2,2); L << 1, 0, 2, 3; > > Vector2d x(1,2), b = L * L.transpose() * x; > > > > LLT<MatrixXd> llt( (MatrixXd()) ); // extra parentheses because C++ ... > > const_cast<MatrixXd&>(llt.matrixLLT())=L; > > > > std::cout << llt.solve(b) << "\n"; > > > > With `LLT<Ref<MatrixXd> >`, the const_cast could be avoided ... > > > > For the LDLT case (and other decompositions), this would require keeping the > > permutations vector at the same size, or strange things will happen. > > And needing const_cast generally is not a sign of a nice API (to say it > > nicely). > > Agreed these solutions sound a bit suboptimal. > > One advantage of leaning on the type system in the TriangularView approach > is it obviates any need for validation of the inputs. We could also avoid > this by simply declaring that the implementation only looks at the > lower/upper triangular portion of the input, and skips validation > altogether. This seems more in the spirit of the current implementation, > which only performs validation in the sense of setting m_info when a > computation goes awry. The situation with permutations in LDLT is more subtle, though.
> If we offer these constructors, we could offer the same > constructors for the standard "copy" variant. What do you mean by "standard "copy" variant"? I've also thought about using triangularView(), but I thought this might be a bit confusing: are you passing the already computed L factor or the triangular-part of the self-adjoint matrix that needs to be factorized? (ok maybe this is a bit far fetched, but I read so many confused user on SO...) Also, I'm not sure it generalizes well to other facotizations. I've also thought about: auto L = A.triangularView<Lower>(); LLT<MatrixXd> llt = L*L.transpose(); but sadly it is not possible to know at compile-time that L is used on both sides, so we would have to check it at runtime; and this is also a bit painful on our side to be sure we catch all cases... Regarding validation of the input, we'll simply document that only one triangular part is read and the rest is ignored. If you want a currently working workaround without const cast: MatrixXd L = ...; MatrixXd tmp = MatrixXd::Identity(L.rows(),L.cols()); LLT<Ref<MatrixXd>> llt(tmp); tmp = L; except that rcond() will be wrong. For LDLT isPositive/isNegative will be wrong too.
(In reply to Gael Guennebaud from comment #7) > What do you mean by "standard "copy" variant"? Essentially, this (and related variants): MatrixXd L; LLT<MatrixXd> llt_wrapper(InitFromFactors, L); i.e., the case without the Ref<> > I've also thought about using triangularView(), but I thought this might be > a bit confusing: are you passing the already computed L factor or the [...] Agreed, especially in generic code, or when we generalize this to other decompositions. Even if there was only one interpretation "that makes sense", it would not be obvious from the API. > Also, I'm not sure it generalizes well to other facotizations. I've also > thought about: > > auto L = A.triangularView<Lower>(); > LLT<MatrixXd> llt = L*L.transpose(); Hm, interesting. The necessary compile-time check would actually very often get optimized away (unless it is really used for X*Y.transpose(), with different X,Y -- but that would also unlikely be symmetric). But it would be bad, if we accidentally miss a valid case and compute the product as well as the decomposition again. Some static method like this could be an alternative to tagged constructors: LLT<MatrixXd> llt = LLT<MatrixXd>::ReconstructFrom(L); Quite a bit longer, but this kind of operation should not be too common. > If you want a currently working workaround without const cast: > > MatrixXd L = ...; > MatrixXd tmp = MatrixXd::Identity(L.rows(),L.cols()); > LLT<Ref<MatrixXd>> llt(tmp); > tmp = L; Yes, you can (and should) even avoid decomposing the identity, by starting with an empty `tmp`. > except that rcond() will be wrong. For LDLT isPositive/isNegative will be > wrong too. I think those will be wrong after rankUpdate(), as well ...
(In reply to Christoph Hertzberg from comment #8) > (In reply to Gael Guennebaud from comment #7) > > auto L = A.triangularView<Lower>(); > > LLT<MatrixXd> llt = L*L.transpose(); > > Hm, interesting. The necessary compile-time check would actually very often > get optimized away (unless it is really used for X*Y.transpose(), with > different X,Y -- but that would also unlikely be symmetric). But it would be > bad, if we accidentally miss a valid case and compute the product as well as > the decomposition again. yes, that's a bit far fetched. > Some static method like this could be an alternative to tagged constructors: > LLT<MatrixXd> llt = LLT<MatrixXd>::ReconstructFrom(L); > Quite a bit longer, but this kind of operation should not be too common. then we also need a copy ctor, and thus operator=, etc. > > If you want a currently working workaround without const cast: > > > > MatrixXd L = ...; > > MatrixXd tmp = MatrixXd::Identity(L.rows(),L.cols()); > > LLT<Ref<MatrixXd>> llt(tmp); > > tmp = L; > > Yes, you can (and should) even avoid decomposing the identity, by starting > with an empty `tmp`. not really because you need to properly set the m_initialized member to true (or always disable assertions :( ), which is why llt.matrixL() = L; (without const_cast) is not sufficient. > > except that rcond() will be wrong. For LDLT isPositive/isNegative will be > > wrong too. > > I think those will be wrong after rankUpdate(), as well ... right, I guess we should track whether those quantities are available or not.
-- GitLab Migration Automatic Message -- This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1599.