add relative_error and project_array functions

docs/source/getting_started.rst

@@ -270,7 +270,7 @@ the H1-product. For this we use the :meth:`~pymor.operators.interfaces.OperatorI
 method:
 
 >>> import numpy as np
->>> gram_matrix = d.h1_product.apply2(reductor.RB, reductor.RB)
+>>> gram_matrix = RB.gramian(d.h1_product)
 >>> print(np.max(np.abs(gram_matrix - np.eye(32))))
 5.86218898944e-14
 
@@ -297,7 +297,7 @@ Finally we compute the reduction error and display the reduced solution along wi
 the detailed solution and the error:
 
 >>> ERR = U - U_red
->>> print(d.h1_norm(ERR))
+>>> print(ERR.norm(d.h1_product))
 [ 0.00944595]
 >>> d.visualize((U, U_red, ERR),
 ...             legend=('Detailed', 'Reduced', 'Error'),

docs/source/technical_overview.rst

@@ -65,7 +65,7 @@ operating on objects of the following types:
     .. |empty|            replace:: :meth:`~pymor.vectorarrays.interfaces.VectorSpaceInterface.empty`
     .. |id|               replace:: :meth:`~pymor.vectorarrays.interfaces.VectorSpaceInterface.id`
     .. |indexed|          replace:: :meth:`indexed <pymor.vectorarrays.interfaces.VectorArrayInterface.__getitem__>`
-    .. |inner products|   replace:: :meth:`inner products <pymor.vectorarrays.interfaces.VectorArrayInterface.dot>`
+    .. |inner products|   replace:: :meth:`inner products <pymor.vectorarrays.interfaces.VectorArrayInterface.inner>`
     .. |lincomb|          replace:: :meth:`~pymor.vectorarrays.interfaces.VectorArrayInterface.lincomb`
     .. |make_array|       replace:: :meth:`~pymor.vectorarrays.interfaces.VectorSpaceInterface.make_array`
     .. |removed|          replace:: :meth:`deleted <pymor.vectorarrays.interfaces.VectorArrayInterface.__delitem__>`

src/pymor/algorithms/basic.py

@@ -67,3 +67,37 @@ def almost_equal(U, V, product=None, norm=None, rtol=1e-14, atol=1e-14):
     ERR_norm = norm(X)
 
     return ERR_norm <= atol + V_norm * rtol
+
+
+def relative_error(U, V, product=None):
+    """Compute error between U and V relative to norm of U."""
+    return (U - V).norm(product) / U.norm(product)
+
+
+def project_array(U, basis, product=None, orthonormal=True):
+    """Orthogonal projection of |VectorArray| onto subspace.
+
+    Parameters
+    ----------
+    U
+        The |VectorArray| to project.
+    basis
+        |VectorArray| of basis vectors for the subspace onto which
+        to project.
+    product
+        Inner product |Operator| w.r.t. which to project.
+    orthonormal
+        If `True`, the vectors in `basis` are assumed to be orthonormal
+        w.r.t. `product`.
+
+    Returns
+    -------
+    The projected |VectorArray|.
+    """
+    if orthonormal:
+        return basis.lincomb(U.inner(basis, product))
+    else:
+        gramian = basis.gramian(product)
+        rhs = basis.inner(U, product)
+        coeffs = np.linalg.solve(gramian, rhs).T
+        return basis.lincomb(coeffs)

src/pymor/algorithms/gram_schmidt.py

@@ -61,10 +61,7 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
     remove = []
     for i in range(offset, len(A)):
         # first calculate norm
-        if product is None:
-            initial_norm = A[i].l2_norm()[0]
-        else:
-            initial_norm = np.sqrt(product.pairwise_apply2(A[i], A[i]))[0]
+        initial_norm = A[i].norm(product)[0]
 
         if initial_norm < atol:
             logger.info("Removing vector {} of norm {}".format(i, initial_norm))
@@ -90,17 +87,11 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
                 for j in range(i):
                     if j in remove:
                         continue
-                    if product is None:
-                        p = A[i].pairwise_dot(A[j])[0]
-                    else:
-                        p = product.pairwise_apply2(A[i], A[j])[0]
+                    p = A[i].pairwise_inner(A[j], product)[0]
                     A[i].axpy(-p, A[j])
 
                 # calculate new norm
-                if product is None:
-                    old_norm, norm = norm, A[i].l2_norm()[0]
-                else:
-                    old_norm, norm = norm, np.sqrt(product.pairwise_apply2(A[i], A[i])[0])
+                old_norm, norm = norm, A[i].norm(product)[0]
 
                 # remove vector if it got too small:
                 if norm / initial_norm < rtol:
@@ -115,10 +106,7 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
         del A[remove]
 
     if check:
-        if product:
-            error_matrix = product.apply2(A[offset:len(A)], A)
-        else:
-            error_matrix = A[offset:len(A)].dot(A)
+        error_matrix = A[offset:len(A)].inner(A, product)
         error_matrix[:len(A) - offset, offset:len(A)] -= np.eye(len(A) - offset)
         if error_matrix.size > 0:
             err = np.max(np.abs(error_matrix))
@@ -176,10 +164,7 @@ def gram_schmidt_biorth(V, W, product=None, reiterate=True, reiteration_threshol
     # main loop
     for i in range(len(V)):
         # calculate norm of V[i]
-        if product is None:
-            initial_norm = V[i].l2_norm()[0]
-        else:
-            initial_norm = np.sqrt(product.pairwise_apply2(V[i], V[i]))[0]
+        initial_norm = V[i].norm(product)[0]
 
         # project V[i]
         if i == 0:
@@ -197,26 +182,17 @@ def gram_schmidt_biorth(V, W, product=None, reiterate=True, reiteration_threshol
 
                 for j in range(i):
                     # project by (I - V[j] * W[j]^T * E)
-                    if product is None:
-                        p = W[j].pairwise_dot(V[i])[0]
-                    else:
-                        p = product.pairwise_apply2(W[j], V[i])[0]
+                    p = W[j].pairwise_inner(V[i], product)[0]
                     V[i].axpy(-p, V[j])
 
                 # calculate new norm
-                if product is None:
-                    old_norm, norm = norm, V[i].l2_norm()[0]
-                else:
-                    old_norm, norm = norm, np.sqrt(product.pairwise_apply2(V[i], V[i])[0])
+                old_norm, norm = norm, V[i].norm(product)[0]
 
             if norm > 0:
                 V[i].scal(1 / norm)
 
         # calculate norm of W[i]
-        if product is None:
-            initial_norm = W[i].l2_norm()[0]
-        else:
-            initial_norm = np.sqrt(product.pairwise_apply2(W[i], W[i]))[0]
+        initial_norm = W[i].norm(product)[0]
 
         # project W[i]
         if i == 0:
@@ -234,33 +210,21 @@ def gram_schmidt_biorth(V, W, product=None, reiterate=True, reiteration_threshol
 
                 for j in range(i):
                     # project by (I - W[j] * V[j]^T * E)
-                    if product is None:
-                        p = V[j].pairwise_dot(W[i])[0]
-                    else:
-                        p = product.pairwise_apply2(V[j], W[i])[0]
+                    p = V[j].pairwise_inner(W[i], product)[0]
                     W[i].axpy(-p, W[j])
 
                 # calculate new norm
-                if product is None:
-                    old_norm, norm = norm, W[i].l2_norm()[0]
-                else:
-                    old_norm, norm = norm, np.sqrt(product.pairwise_apply2(W[i], W[i])[0])
+                old_norm, norm = norm, W[i].norm(product)[0]
 
             if norm > 0:
                 W[i].scal(1 / norm)
 
         # rescale V[i]
-        if product is None:
-            p = W[i].pairwise_dot(V[i])[0]
-        else:
-            p = product.pairwise_apply2(W[i], V[i])[0]
+        p = W[i].pairwise_inner(V[i], product)[0]
         V[i].scal(1 / p)
 
     if check:
-        if product:
-            error_matrix = product.apply2(W, V)
-        else:
-            error_matrix = W.dot(V)
+        error_matrix = W.inner(V, product)
         error_matrix -= np.eye(len(V))
         if error_matrix.size > 0:
             err = np.max(np.abs(error_matrix))

src/pymor/algorithms/pod.py

@@ -72,7 +72,7 @@ def pod(A, modes=None, product=None, rtol=4e-8, atol=0., l2_err=0.,
     logger = getLogger('pymor.algorithms.pod.pod')
 
     with logger.block('Computing Gramian ({} vectors) ...'.format(len(A))):
-        B = A.gramian() if product is None else product.apply2(A, A)
+        B = A.gramian(product)
 
         if symmetrize:     # according to rbmatlab this is necessary due to rounding
             B = B + B.T
@@ -111,13 +111,8 @@ def pod(A, modes=None, product=None, rtol=4e-8, atol=0., l2_err=0.,
 
     if check:
         logger.info('Checking orthonormality ...')
-        if not product and not float_cmp_all(POD.dot(POD), np.eye(len(POD)),
-                                             atol=check_tol, rtol=0.):
-            err = np.max(np.abs(POD.dot(POD) - np.eye(len(POD))))
-            raise AccuracyError('result not orthogonal (max err={})'.format(err))
-        elif product and not float_cmp_all(product.apply2(POD, POD), np.eye(len(POD)),
-                                           atol=check_tol, rtol=0.):
-            err = np.max(np.abs(product.apply2(POD, POD) - np.eye(len(POD))))
+        if not float_cmp_all(POD.inner(POD, product), np.eye(len(POD)), atol=check_tol, rtol=0.):
+            err = np.max(np.abs(POD.inner(POD, product) - np.eye(len(POD))))
             raise AccuracyError('result not orthogonal (max err={})'.format(err))
         if len(POD) < len(EVECS):
             raise AccuracyError('additional orthonormalization removed basis vectors')

src/pymor/algorithms/projection.py

@@ -87,10 +87,7 @@ def action_ZeroOperator(self, op, range_basis, source_basis, product=None):
     @match_class(ConstantOperator)
     def action_ConstantOperator(self, op, range_basis, source_basis, product=None):
         if range_basis is not None:
-            if product:
-                projected_value = NumpyVectorSpace.make_array(product.apply2(range_basis, op._value).T, op.range.id)
-            else:
-                projected_value = NumpyVectorSpace.make_array(range_basis.dot(op._value).T, op.range.id)
+            projected_value = NumpyVectorSpace.make_array(range_basis.inner(op._value, product).T, op.range.id)
         else:
             projected_value = op._value
         if source_basis is None:
@@ -184,13 +181,9 @@ def action_EmpiricalInterpolatedOperator(self, op, range_basis, source_basis, pr
             raise RuleNotMatchingError('Has no restricted operator or source_basis is None')
         else:
             if range_basis is not None:
-                if product is None:
-                    projected_collateral_basis = NumpyVectorSpace.make_array(op.collateral_basis.dot(range_basis),
-                                                                             op.range.id)
-                else:
-                    projected_collateral_basis = NumpyVectorSpace.make_array(product.apply2(op.collateral_basis,
-                                                                                            range_basis),
-                                                                             op.range.id)
+                projected_collateral_basis = NumpyVectorSpace.make_array(op.collateral_basis.inner(range_basis,
+                                                                                                   product),
+                                                                         op.range.id)
             else:
                 projected_collateral_basis = op.collateral_basis
 

src/pymor/basic.py

@@ -8,7 +8,7 @@
 to have the most important parts of pyMOR directly available.
 """
 
-from pymor.algorithms.basic import almost_equal
+from pymor.algorithms.basic import almost_equal, relative_error, project_array
 from pymor.algorithms.ei import interpolate_operators, ei_greedy, deim
 from pymor.algorithms.error import reduction_error_analysis
 from pymor.algorithms.gram_schmidt import gram_schmidt, gram_schmidt_biorth

src/pymor/vectorarrays/interfaces.py

@@ -230,6 +230,17 @@ def dot(self, other):
         """
         pass
 
+    def inner(self, other, product=None):
+        """Inner products w.r.t. given product |Operator|.
+
+        Equivalent to `self.dot(other)` if `product` is None,
+        else equivalent to `product.apply2(self, other)`.
+        """
+        if product is None:
+            return self.dot(other)
+        else:
+            return product.apply2(self, other)
+
     @abstractmethod
     def pairwise_dot(self, other):
         """Returns the pairwise inner products between |VectorArray| elements.
@@ -248,6 +259,17 @@ def pairwise_dot(self, other):
         """
         pass
 
+    def pairwise_inner(self, other, product=None):
+        """Pairwise inner products w.r.t. given product |Operator|.
+
+        Equivalent to `self.pairwise_dot(other)` if `product` is None,
+        else equivalent to `product.pairwise_apply2(self, other)`.
+        """
+        if product is None:
+            return self.pairwise_dot(other)
+        else:
+            return product.pairwise_apply2(self, other)
+
     @abstractmethod
     def lincomb(self, coefficients):
         """Returns linear combinations of the vectors contained in the array.
@@ -272,6 +294,28 @@ def lincomb(self, coefficients):
         """
         pass
 
+    def norm(self, product=None):
+        """Norm w.r.t. given inner product |Operator|.
+
+        Equivalent to `self.l2_norm()` if `product` is None,
+        else equivalent to `np.sqrt(product.pairwise_apply2(self, self))`.
+        """
+        if product is None:
+            return self.l2_norm()
+        else:
+            return np.sqrt(product.pairwise_apply2(self, self))
+
+    def norm2(self, product=None):
+        """Squared norm w.r.t. given inner product |Operator|.
+
+        Equivalent to `self.l2_norm2()` if `product` is None,
+        else equivalent to `product.pairwise_apply2(self, self)`.
+        """
+        if product is None:
+            return self.l2_norm2()
+        else:
+            return product.pairwise_apply2(self, self)
+
     @abstractmethod
     def l1_norm(self):
         """The l1-norms of the vectors contained in the array.
@@ -351,9 +395,9 @@ def amax(self):
         """
         pass
 
-    def gramian(self):
-        """Shorthand for `self.dot(self)`."""
-        return self.dot(self)
+    def gramian(self, product=None):
+        """Shorthand for `self.inner(self, product)`."""
+        return self.inner(self, product)
 
     def __add__(self, other):
         """The pairwise sum of two |VectorArrays|."""

src/pymor/vectorarrays/list.py

@@ -320,7 +320,9 @@ def pairwise_dot(self, other):
         assert len(self._list) == len(other)
         return np.array([a.dot(b) for a, b in zip(self._list, other._list)])
 
-    def gramian(self):
+    def gramian(self, product=None):
+        if product is not None:
+            return super().gramian(product)
         l = len(self._list)
         R = np.empty((l, l))
         for i in range(l):

src/pymortests/algorithms/basic.py

@@ -7,8 +7,11 @@
 import pytest
 import numpy as np
 
-from pymor.algorithms.basic import almost_equal
+from pymor.algorithms.basic import almost_equal, project_array, relative_error
+from pymor.algorithms.gram_schmidt import gram_schmidt
 from pymor.operators.constructions import induced_norm
+from pymor.operators.numpy import NumpyMatrixOperator
+from pymor.vectorarrays.numpy import NumpyVectorSpace
 from pymortests.fixtures.vectorarray import \
     (vector_array_without_reserve, vector_array, compatible_vector_array_pair_without_reserve,
      compatible_vector_array_pair, incompatible_vector_array_pair)
@@ -169,3 +172,25 @@ def test_almost_equal_wrong_ind(compatible_vector_array_pair):
             c1, c2 = v1.copy(), v2.copy()
             with pytest.raises(Exception):
                 almost_equal(c1[ind], c2[ind], norm=n)
+
+
+def test_project_array():
+    np.random.seed(123)
+    U = NumpyVectorSpace.from_data(np.random.random((2, 10)))
+    basis = NumpyVectorSpace.from_data(np.random.random((3, 10)))
+    U_p = project_array(U, basis, orthonormal=False)
+    onb = gram_schmidt(basis)
+    U_p2 = project_array(U, onb, orthonormal=True)
+    assert np.all(relative_error(U_p, U_p2) < 1e-10)
+
+
+def test_project_array_with_product():
+    np.random.seed(123)
+    U = NumpyVectorSpace.from_data(np.random.random((1, 10)))
+    basis = NumpyVectorSpace.from_data(np.random.random((3, 10)))
+    product = np.random.random((10, 10))
+    product = NumpyMatrixOperator(product.T.dot(product))
+    U_p = project_array(U, basis, product=product, orthonormal=False)
+    onb = gram_schmidt(basis, product=product)
+    U_p2 = project_array(U, onb, product=product, orthonormal=True)
+    assert np.all(relative_error(U_p, U_p2, product) < 1e-10)