Commit 7b482693 authored by Stephan Rave's avatar Stephan Rave Committed by GitHub

Merge pull request #885 from pymor/vectorarrayoperator_apply_inverse

Implement VectorArrayOperator.apply_inverse
parents 4793ad25 35e47852
Pipeline #51142 passed with stages
in 55 minutes and 22 seconds
......@@ -841,6 +841,24 @@ class VectorArrayOperator(Operator):
else:
return self.range.make_array(self.array.dot(U).T)
def apply_inverse(self, V, mu=None, least_squares=False):
if not least_squares and len(self.array) != self.array.dim:
raise InversionError
from pymor.algorithms.gram_schmidt import gram_schmidt
from numpy.linalg import lstsq
Q, R = gram_schmidt(self.array, return_R=True, reiterate=False)
if self.adjoint:
v = lstsq(R.T.conj(), V.to_numpy().T)[0]
U = Q.lincomb(v.T)
else:
v = Q.dot(V)
u = lstsq(R, v)[0]
U = self.source.make_array(u.T)
return U
def apply_adjoint(self, V, mu=None):
assert V in self.range
if not self.adjoint:
......
......@@ -11,7 +11,7 @@ from pymor.algorithms.to_matrix import to_matrix
from pymor.core.exceptions import InversionError, LinAlgError
from pymor.operators.block import BlockDiagonalOperator
from pymor.operators.constructions import (SelectionOperator, InverseOperator, InverseAdjointOperator, IdentityOperator,
LincombOperator)
LincombOperator, VectorArrayOperator)
from pymor.operators.numpy import NumpyMatrixOperator
from pymor.parameters.base import ParameterType
from pymor.parameters.functionals import GenericParameterFunctional, ExpressionParameterFunctional
......@@ -417,3 +417,36 @@ def test_InverseAdjointOperator(operator_with_arrays):
rtol=rtol, atol=atol))
except (InversionError, LinAlgError, NotImplementedError):
pass
def test_vectorarray_op_apply_inverse():
np.random.seed(1234)
O = np.random.random((5, 5))
op = VectorArrayOperator(NumpyVectorSpace.make_array(O))
V = op.range.random()
U = op.apply_inverse(V)
v = V.to_numpy()
u = np.linalg.solve(O.T, v.ravel())
assert np.all(almost_equal(U, U.space.from_numpy(u)))
def test_vectorarray_op_apply_inverse_lstsq():
np.random.seed(1234)
O = np.random.random((3, 5))
op = VectorArrayOperator(NumpyVectorSpace.make_array(O))
V = op.range.random()
U = op.apply_inverse(V, least_squares=True)
v = V.to_numpy()
u = np.linalg.lstsq(O.T, v.ravel())[0]
assert np.all(almost_equal(U, U.space.from_numpy(u)))
def test_adjoint_vectorarray_op_apply_inverse_lstsq():
np.random.seed(1234)
O = np.random.random((3, 5))
op = VectorArrayOperator(NumpyVectorSpace.make_array(O), adjoint=True)
V = op.range.random()
U = op.apply_inverse(V, least_squares=True)
v = V.to_numpy()
u = np.linalg.lstsq(O, v.ravel())[0]
assert np.all(almost_equal(U, U.space.from_numpy(u)))
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment