# Source code for pymor.algorithms.gram_schmidt

# This file is part of the pyMOR project (http://www.pymor.org).

import numpy as np

from pymor.core.defaults import defaults
from pymor.core.exceptions import AccuracyError
from pymor.core.logger import getLogger

[docs]@defaults('atol', 'rtol', 'reiterate', 'reiteration_threshold', 'check', 'check_tol')
def gram_schmidt(A, product=None, return_R=False, atol=1e-13, rtol=1e-13, offset=0,
reiterate=True, reiteration_threshold=1e-1, check=True, check_tol=1e-3,
copy=True):
"""Orthonormalize a |VectorArray| using the modified Gram-Schmidt algorithm.

Parameters
----------
A
The |VectorArray| which is to be orthonormalized.
product
The inner product |Operator| w.r.t. which to orthonormalize.
If None, the Euclidean product is used.
return_R
If True, the R matrix from QR decomposition is returned.
atol
Vectors of norm smaller than atol are removed from the array.
rtol
Relative tolerance used to detect linear dependent vectors
(which are then removed from the array).
offset
Assume that the first offset vectors are already orthonormal and start the
algorithm at the offset + 1-th vector.
reiterate
If True, orthonormalize again if the norm of the orthogonalized vector is
much smaller than the norm of the original vector.
reiteration_threshold
If reiterate is True, re-orthonormalize if the ratio between the norms of
the orthogonalized vector and the original vector is smaller than this value.
check
If True, check if the resulting |VectorArray| is really orthonormal.
check_tol
Tolerance for the check.
copy
If True, create a copy of A instead of modifying A in-place.

Returns
-------
Q
The orthonormalized |VectorArray|.
R
The upper-triangular/trapezoidal matrix (if compute_R is True).
"""

logger = getLogger('pymor.algorithms.gram_schmidt.gram_schmidt')

if copy:
A = A.copy()

# main loop
R = np.eye(len(A))
remove = []  # indices of to be removed vectors
for i in range(offset, len(A)):
# first calculate norm
initial_norm = A[i].norm(product)

if initial_norm < atol:
logger.info(f"Removing vector {i} of norm {initial_norm}")
remove.append(i)
continue

if i == 0:
A.scal(1 / initial_norm)
R[i, i] = initial_norm
else:
norm = initial_norm
# If reiterate is True, reiterate as long as the norm of the vector changes
# strongly during orthogonalization (due to Andreas Buhr).
while True:
# orthogonalize to all vectors left
for j in range(i):
if j in remove:
continue
p = A[j].pairwise_inner(A[i], product)
A[i].axpy(-p, A[j])
common_dtype = np.promote_types(R.dtype, type(p))
R = R.astype(common_dtype, copy=False)
R[j, i] += p

# calculate new norm
old_norm, norm = norm, A[i].norm(product)

# remove vector if it got too small
if norm < rtol * initial_norm:
logger.info(f"Removing linearly dependent vector {i}")
remove.append(i)
break

# check if reorthogonalization should be done
if reiterate and norm < reiteration_threshold * old_norm:
logger.info(f"Orthonormalizing vector {i} again")
else:
A[i].scal(1 / norm)
R[i, i] = norm
break

if remove:
del A[remove]
R = np.delete(R, remove, axis=0)

if check:
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))
if err >= check_tol:
raise AccuracyError(f"result not orthogonal (max err={err})")

if return_R:
return A, R
else:
return A

[docs]def gram_schmidt_biorth(V, W, product=None,
reiterate=True, reiteration_threshold=1e-1, check=True, check_tol=1e-3,
copy=True):
"""Biorthonormalize a pair of |VectorArrays| using the biorthonormal Gram-Schmidt process.

See Algorithm 1 in [BKS11]_.

Parameters
----------
V, W
The |VectorArrays| which are to be biorthonormalized.
product
The inner product |Operator| w.r.t. which to biorthonormalize.
If None, the Euclidean product is used.
reiterate
If True, orthonormalize again if the norm of the orthogonalized vector is
much smaller than the norm of the original vector.
reiteration_threshold
If reiterate is True, re-orthonormalize if the ratio between the norms of
the orthogonalized vector and the original vector is smaller than this value.
check
If True, check if the resulting |VectorArray| is really orthonormal.
check_tol
Tolerance for the check.
copy
If True, create a copy of V and W instead of modifying V and W in-place.

Returns
-------
The biorthonormalized |VectorArrays|.
"""
assert V.space == W.space
assert len(V) == len(W)

logger = getLogger('pymor.algorithms.gram_schmidt.gram_schmidt_biorth')

if copy:
V = V.copy()
W = W.copy()

# main loop
for i in range(len(V)):
# calculate norm of V[i]
initial_norm = V[i].norm(product)

# project V[i]
if i == 0:
V.scal(1 / initial_norm)
else:
norm = initial_norm
# If reiterate is True, reiterate as long as the norm of the vector changes
# strongly during projection.
while True:
for j in range(i):
# project by (I - V[j] * W[j]^T * E)
p = W[j].pairwise_inner(V[i], product)
V[i].axpy(-p, V[j])

# calculate new norm
old_norm, norm = norm, V[i].norm(product)

# check if reorthogonalization should be done
if reiterate and norm < reiteration_threshold * old_norm:
logger.info(f"Projecting vector V[{i}] again")
else:
V[i].scal(1 / norm)
break

# calculate norm of W[i]
initial_norm = W[i].norm(product)

# project W[i]
if i == 0:
W.scal(1 / initial_norm)
else:
norm = initial_norm
# If reiterate is True, reiterate as long as the norm of the vector changes
# strongly during projection.
while True:
for j in range(i):
# project by (I - W[j] * V[j]^T * E)
p = V[j].pairwise_inner(W[i], product)
W[i].axpy(-p, W[j])

# calculate new norm
old_norm, norm = norm, W[i].norm(product)

# check if reorthogonalization should be done
if reiterate and norm < reiteration_threshold * old_norm:
logger.info(f"Projecting vector W[{i}] again")
else:
W[i].scal(1 / norm)
break

# rescale V[i]
p = W[i].pairwise_inner(V[i], product)
V[i].scal(1 / p)

if check:
error_matrix = W.inner(V, product)
error_matrix -= np.eye(len(V))
if error_matrix.size > 0:
err = np.max(np.abs(error_matrix))
if err >= check_tol:
raise AccuracyError(f"result not biorthogonal (max err={err})")

return V, W