Cholesky decompositions of `LinearOperator`s #618

marvinpfoertner · 2022-01-23T19:28:57Z

In a Nutshell

This PR implements Cholesky decompositions for ProbNum's linops library.
To achieve this goal, I also had to implement symmetrization and parts of a matrix property system.

Detailed Description

LinearOperator.cholesky(lower: bool) -> LinearOperator
LinearOperator.symmetrize() -> LinearOperator implementing "1/2 * (A + A^T)"
Matrix properties
- is_symmetric
- is_{upper,lower}_triangular
- is_positive_definite
Fixed bugs in Kronecker.{inv,cond}
_InverseLinearOperator uses Cholesky factorization for SPD matrices
Fixed bug in _InverseLinearOperator of singular matrix
Added new tests and test cases

Related Issues

Implements parts of #569 and #571

codecov · 2022-01-23T19:46:21Z

Codecov Report

Merging #618 (285b304) into main (846f0bc) will increase coverage by 0.11%.
The diff coverage is 90.90%.

@@            Coverage Diff             @@
##             main     #618      +/-   ##
==========================================
+ Coverage   89.40%   89.52%   +0.11%     
==========================================
  Files         181      181              
  Lines        6684     6875     +191     
  Branches     1026     1064      +38     
==========================================
+ Hits         5976     6155     +179     
- Misses        479      488       +9     
- Partials      229      232       +3

Impacted Files	Coverage Δ
src/probnum/linops/_linear_operator.py	`86.15% <87.50%> (+0.67%)`	⬆️
src/probnum/linops/_scaling.py	`82.55% <89.47%> (+0.85%)`	⬆️
src/probnum/linops/_arithmetic_fallbacks.py	`94.44% <100.00%> (+1.58%)`	⬆️
src/probnum/linops/_kronecker.py	`90.94% <100.00%> (+3.32%)`	⬆️

JonathanWenger

It's great to finally have this property system. I have a few comments which should be easily resolved. Thanks for the effort! One question to think about for the future (not this PR) what do we do with semi-definite operators? Add another property _is_positive_semi_definite? Often there are algorithms which also work if the linear operator is at least positive semi-definite.

src/probnum/linops/_linear_operator.py

JonathanWenger · 2022-01-28T17:33:57Z

src/probnum/linops/_linear_operator.py

+        numpy.linalg.LinAlgError
+            If the :class:`LinearOperator` is not positive definite.
+        """
+        if not self.is_symmetric:


What happens if the is_symmetric and is_positive_definite properties are None, i.e. it is unknown whether the matrix is spd? According to the docstring .cholesky() should still run and if it succeeds, set the properties accordingly. However if self.is_symmetric is None, this conditional should throw an exception, if i'm not mistaken.

According to the Raises section of the docstring, the method raises an exception if is_symmetric is not True, i.e. the method will fail if it is not known that the operator is indeed symmetric.
Similarly I state that the operator raises an exception if the matrix is not positive definite, while still setting linop.is_positive_definite = False.
This is intentional as it can be used as a test for positive definiteness.
So I don't really see a contradiction to the docstring here. Or am I missing something?

My mistake! We might want to be lenient here and only raise an exception if the properties are False. Otherwise this may lead to a lot of boilerplate in particular when prototyping. Your call.

In fact, I think it would be more in line with Python philosophy to attempt the Cholesky decomposition if it is not known whether the linop is spd and then raise an exception should Cholesky fail.

In principle, I agree with the second point. In fact this is what happens if is_positive_definite is None: We attempt the Cholesky and if it fails, we set it to False and raise an exception. However, after giving this some thought, this issue is much more complicated with is_symmetric.
Since matrices are rarely exactly symmetric, an automatic check of symmetry would be much too restrictive, unless we have some sort of "asymmetry tolerance", which is just another one of these arbitrary parameters to make things work. Personally, I came to the conclusion that it is better if the user has to call .symmetrize on the operator before calling .cholesky because it makes this whole symmetry business much more transparent.
Moreover, if symmetrize is too costly, the user can always choose to set linop.symmetric = True manually, even if the property is only fulfilled approximately.

Note that NumPy and SciPy circumvent this issue by only operating on the lower triangle of the matrix.

>>> a = np.array([ ... [1.0, 2.0], ... [0.9, 1.0], ... ]) >>> np.linalg.cholesky(a) array([[1. , 0. ], [0.9 , 0.43588989]]) >>> scipy.linalg.cholesky(a, lower=True) array([[1. , 0. ], [0.9 , 0.43588989]])

IMHO, this is quite dangerous, since it can produce very subtle bugs.

Sounds good to me! Thanks for the detailed explanation.

If we find a good solution, I'm also fine with removing this in the future.
However, I would propose to make this a separate PR, since this would increase complexity by quite a bit.
We can always make interfaces less restrictive

JonathanWenger · 2022-01-28T17:42:07Z

src/probnum/linops/_linear_operator.py

@@ -520,6 +709,59 @@ def inv(self) -> "LinearOperator":

        return self.__inverse()

+    def symmetrize(self) -> LinearOperator:


By introduction of this method the linear operator Symmetrize in _kronecker.py should just be removed, I think. Otherwise there are two different ways of achieving the same thing.

Symmetrize in _kronecker.py symmetrizes vectorized matrices used in Kronecker algebra, so it operates "on the arguments" of linear operators, while linop.symmetrize symmetrizes the LinearOperator itself, so I'd say they serve distinct purposes.

They do, but Symmetrize is really only used for SymmetricKronecker and is arguably confusing in its naming if there is a symmetrize() method. My suggestion would be to fold Symmetrize into the SymmetricKronecker implementation as a (future) refactor. I'm okay with considering this as out of scope for this PR.

I agree, let's make a separate PR here 👍

src/probnum/linops/_scaling.py

marvinpfoertner · 2022-01-28T18:45:29Z

It's great to finally have this property system. I have a few comments which should be easily resolved. Thanks for the effort! One question to think about for the future (not this PR) what do we do with semi-definite operators? Add another property _is_positive_semi_definite? Often there are algorithms which also work if the linear operator is at least positive semi-definite.

I agree, we should add more matrix properties over time. I also think that we will need is_indefinite, as well as is_diagonal and is_orthogonal (the latter for SVDs and eigendecompositions). I started a list of these properties in #571 and #569. Feel free to add is_positive_semidefinite there.

JonathanWenger · 2022-01-28T18:46:24Z

src/probnum/linops/_linear_operator.py

+        numpy.linalg.LinAlgError
+            If the :class:`LinearOperator` is not positive definite.
+        """
+        if not self.is_symmetric:


My mistake! We might want to be lenient here and only raise an exception if the properties are False. Otherwise this may lead to a lot of boilerplate in particular when prototyping. Your call.

JonathanWenger · 2022-01-28T18:54:43Z

src/probnum/linops/_linear_operator.py

@@ -520,6 +709,59 @@ def inv(self) -> "LinearOperator":

        return self.__inverse()

+    def symmetrize(self) -> LinearOperator:


They do, but Symmetrize is really only used for SymmetricKronecker and is arguably confusing in its naming if there is a symmetrize() method. My suggestion would be to fold Symmetrize into the SymmetricKronecker implementation as a (future) refactor. I'm okay with considering this as out of scope for this PR.

marvinpfoertner added improvement Improvements of existing functionality linops Issues related to linear operators labels Jan 23, 2022

marvinpfoertner self-assigned this Jan 23, 2022

marvinpfoertner force-pushed the linop-cholesky branch 2 times, most recently from 6d65351 to a50d087 Compare January 24, 2022 13:03

marvinpfoertner changed the title ~~[Draft] Cholesky decompositions of LinearOperators~~ Cholesky decompositions of LinearOperators Jan 24, 2022

marvinpfoertner marked this pull request as ready for review January 24, 2022 17:37

marvinpfoertner requested a review from JonathanWenger as a code owner January 24, 2022 17:37

JonathanWenger requested changes Jan 28, 2022

View reviewed changes

marvinpfoertner force-pushed the linop-cholesky branch from 49fc7d2 to 66dfdb6 Compare January 28, 2022 18:50

marvinpfoertner added 16 commits January 28, 2022 19:51

Matrix properties for the linear operator

692dd39

Linop cholesky decomposition

cc50ec1

Kronecker cholesky

a871230

Symmetrization

16a05c4

Add scaling and arithmetic fallback properties

ccb2cd6

Scaling Cholesky

46ee01d

Identity Kronecker Cholesky

dfa74b6

use cholesky for solving

55c7b3f

Inverse linop properties

5104933

Cholesky implies positive definite

7bc0077

Add tags to Scaling test cases

19c02aa

Bugfix for matrix properties

e8d58be

Add simple Cholesky test

8c00be1

Bugfix in Scaling

db6c6ad

Add square tag in Scaling

f8e1ce3

Kronecker SPD cases

1aa6a4d

marvinpfoertner added 25 commits January 28, 2022 19:51

Add test cases for linop fallback arithmetic

92a4a4e

Bugfix in arithmetic fallbacks

591b3f3

More arithmetic cases

76c267a

Asymmetric IdentityKronecker test case

cd7a5cf

Bugfix in Kronecker

e96f63b

New Kronecker test case

eca7e6d

Bugfix in InverseLinearOperator

e2e6a6a

Bugfixes in Kronecker

db3a1f7

Fixes in linop tests

f07f4c6

pylint fixes

14b2530

Test case fixes

b76aff6

Kronecker test fix

75d9f6e

Improve Kronecker cases

7118127

Better Kronecker test case

50e7969

More Cholesky tests

f7c6f04

Improve scaling cases

12454cf

Docstrings and consistency checks for matrix properties

b02dafa

Own testing module for decompositions

d9fed1d

LinearOperator.cholesky documentation

0a0b194

Symmetrize documentation

f587709

Bugfix in the tests

6aa80b7

PyLint fix

957e14a

Updates to the documentation

f3c1a4a

docstring fixes

874b9ac

Incorporate review comments

285b304

marvinpfoertner force-pushed the linop-cholesky branch from 66dfdb6 to 285b304 Compare January 28, 2022 18:51

JonathanWenger approved these changes Jan 28, 2022

View reviewed changes

marvinpfoertner merged commit 6659bb6 into probabilistic-numerics:main Jan 28, 2022

marvinpfoertner deleted the linop-cholesky branch January 28, 2022 19:29

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Cholesky decompositions of `LinearOperator`s #618

Cholesky decompositions of `LinearOperator`s #618

marvinpfoertner commented Jan 23, 2022 •

edited

Loading

codecov bot commented Jan 23, 2022 •

edited

Loading

JonathanWenger left a comment

JonathanWenger Jan 28, 2022

marvinpfoertner Jan 28, 2022

JonathanWenger Jan 28, 2022

JonathanWenger Jan 28, 2022

marvinpfoertner Jan 28, 2022

marvinpfoertner Jan 28, 2022

JonathanWenger Jan 28, 2022

marvinpfoertner Jan 28, 2022

JonathanWenger Jan 28, 2022

marvinpfoertner Jan 28, 2022

JonathanWenger Jan 28, 2022

marvinpfoertner Jan 28, 2022 •

edited

Loading

marvinpfoertner commented Jan 28, 2022

JonathanWenger Jan 28, 2022

JonathanWenger Jan 28, 2022

		@@ -520,6 +709,59 @@ def inv(self) -> "LinearOperator":

		return self.__inverse()

		def symmetrize(self) -> LinearOperator:

Cholesky decompositions of LinearOperators #618

Cholesky decompositions of LinearOperators #618

Conversation

marvinpfoertner commented Jan 23, 2022 • edited Loading

In a Nutshell

Detailed Description

Related Issues

codecov bot commented Jan 23, 2022 • edited Loading

Codecov Report

JonathanWenger left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

marvinpfoertner Jan 28, 2022 • edited Loading

Choose a reason for hiding this comment

marvinpfoertner commented Jan 28, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Cholesky decompositions of `LinearOperator`s #618

Cholesky decompositions of `LinearOperator`s #618

marvinpfoertner commented Jan 23, 2022 •

edited

Loading

codecov bot commented Jan 23, 2022 •

edited

Loading

marvinpfoertner Jan 28, 2022 •

edited

Loading