Caching of the LU decomposition computed by _InverseLinearOperator

[marvinpfoertner]

In a Nutshell

Currently, executing

A = pn.linops.Matrix(...)
x1 = A.inv() @ b1
x2 = A.inv() @ b2

runs scipy.linalg.lu_factor on A.dense twice.
In principle, it would be possible to cache the return value of A.inv(), but that would introduce cyclic references,
which might cause memory overhead before being collected by the garbage collector.
An elegant solution is to cache all expensive quantities in the linear operator A itself, and then create "throwaway" objects like _InverseLinearOperator and TransposeLinearOperator, which are essentially views of A, accessing these quantities through the base linear operator.

Detailed Description

• caching of LU decomposition in LinearOperators
• bugfix in _InverseLinearOperator.__rmatmul__
• add test case for _InverseLinearOperator
• Replace custom implementation of LinearOperator.broadcast_matmul with np.vectorize-based implementation

Related Issues

None

[marvinpfoertner]

@JonathanWenger This PR should be reviewed after #735 is merged, since it branches off of it.

[codecov]

Codecov Report

Merging #737 (f206691) into main (ce4bbba) will decrease coverage by 0.04%.
The diff coverage is 79.16%.

@@            Coverage Diff             @@
##             main     #737      +/-   ##
==========================================
- Coverage   90.40%   90.36%   -0.05%     
==========================================
  Files         197      197              
  Lines        7454     7431      -23     
  Branches      966      964       -2     
==========================================
- Hits         6739     6715      -24     
- Misses        481      482       +1     
  Partials      234      234              

Impacted Files	Coverage Δ
src/probnum/linops/_linear_operator.py	87.38% <79.16%> (-0.67%)	⬇️

[JonathanWenger]

Nice! Elegant solution of the circular dependency issue.

[JonathanWenger]

Suggested change

	The solutions :math:`A^{-1} X` of the linear systems.
	The solution(s) :math:`A^{-1} X` of the linear system.

[marvinpfoertner]

We actually always have multiple solutions here, since this method only works if x.ndim == 2.
I added an appropriate assertion.