removing redundant code from module varprodmd.py

Author: greinerth

pydmd/varprodmd.py

@@ -17,15 +17,17 @@
 """
 
 import warnings
-from typing import Any, Dict, Tuple, Union, Optional
 from types import MappingProxyType
+from typing import Any, Dict, Tuple, Union
+
 import numpy as np
-from scipy.optimize import OptimizeResult, least_squares
 from scipy.linalg import qr
+from scipy.optimize import OptimizeResult, least_squares
+
 from .dmd import DMDBase
 from .dmdoperator import DMDOperator
 from .snapshots import Snapshots
-
+from .utils import compute_svd
 
 OPT_DEF_ARGS = MappingProxyType(
     {
@@ -40,110 +42,6 @@
 )
 
 
-def _svht(
-    sigma_svd: np.ndarray,
-    rows: int,
-    cols: int,
-    sigma: Optional[float] = None,
-) -> int:
-    """
-    Determine optimal rank for svd matrix approximation,
-    based on https://arxiv.org/pdf/1305.5870.pdf.
-
-    :param sigma_svd: Diagonal matrix of "enonomy" SVD as 1d array.
-    :type sigma_svd: np.ndarray
-    :param rows: Number of rows of original data matrix.
-    :type cols: int
-    :param cols: Number of columns of original data matrix.
-    :type cols: int
-    :param sigma: Signal noise level if known, defaults to None.
-    :type sigma: Optional[float], optional
-    :raises ValueError: sigma_svd must be a 1d-array else exception is raised.
-    :return: Optimal rank.
-    :rtype: int
-    """
-
-    if len(sigma_svd.shape) != 1:
-        raise ValueError("Expected 1d array for diagonal matrix!")
-
-    beta = (
-        float(cols) / float(rows) if rows > cols else float(rows) / float(cols)
-    )
-
-    tau_star = 0
-
-    if sigma is not None:
-        sigma = abs(sigma)
-        lambda_star = np.sqrt(
-            2 * (beta + 1)
-            + (8 * beta / (beta + 1 + np.sqrt(beta * beta + 14 * beta + 1)))
-        )
-        tau_star = (
-            lambda_star
-            * sigma
-            * np.sqrt(float(cols) if cols > rows else float(rows))
-        )
-    else:
-        median = np.median(sigma_svd)
-        beta_square = beta * beta
-        beta_cubic = beta * beta_square
-        omega = 0.56 * beta_cubic - 0.95 * beta_square + 1.82 * beta + 1.43
-        tau_star = median * omega
-
-    rank = np.where(sigma_svd >= tau_star)[0]
-
-    if rank.size == 0:
-        return sigma_svd.shape[-1]
-
-    return rank[-1] + 1
-
-
-def _compute_rank(
-    sigma_x: np.ndarray,
-    rows: int,
-    cols: int,
-    rank: Union[float, int],
-    sigma: float = None,
-) -> int:
-    r"""
-    Compute rank without duplicate SVD computation.
-
-    :param sigma_x: Diagonal matrix of SVD as 1d array.
-    :type sigma_x: np.ndarray
-    :param rows: Number of rows of the original data matrix.
-    :type rows: int
-    :param cols: Number of columns of the original data matrix.
-    :type cols: int
-    :param rank: Desired input rank :math:`r`. If rank is an integer,
-        and :math:`r > 0`, the desired rank is used iff possible.
-        This depends on the size of parameter sigma_x. If the desired
-        rank exceeds the size of sigma_x, then the minimum of desired rank
-        and the size of sigma_x is used. If the rank is a float
-        and :math:`0 < r < 1`, than the rank is determined by investigating
-        the cumulative energy of the singular values.
-    :type rank: Union[float, int]
-    :param sigma: Noise level for SVHT algorithm. If it is unknown
-        keep it as None. Defaults to None.
-    :type sigma: float, optional
-    :raises ValueError: ValueError is raised if rank is negative
-        (:math:`r < 0`).
-    :return: Optimal rank.
-    :rtype: int
-    """
-
-    if 0 < rank < 1:
-        sigma_x_square = np.square(sigma_x)
-        cumulative_energy = np.cumsum(sigma_x_square / sigma_x_square.sum())
-        __rank = np.searchsorted(cumulative_energy, rank) + 1
-    elif rank == 0:
-        __rank = _svht(sigma_x, rows, cols, sigma)
-    elif rank >= 1:
-        __rank = int(rank)
-    else:
-        raise ValueError(f"Invalid rank specified, provided {rank}!")
-    return min(__rank, sigma_x.size)
-
-
 def _compute_dmd_ev(
     x_current: np.ndarray,
     x_next: np.ndarray,
@@ -167,22 +65,13 @@ def _compute_dmd_ev(
     :rtype: np.ndarray
     """
 
-    u_x, sigma_x, v_x_t = np.linalg.svd(
-        x_current, full_matrices=False, hermitian=False
-    )
-    __rank = _compute_rank(
-        sigma_x, x_current.shape[0], x_current.shape[1], rank
-    )
+    u_x, sigma_x, v_x = compute_svd(x_current, rank)
 
     # columns of v need to be multiplicated with inverse sigma
-    sigma_inv_approx = np.reciprocal(sigma_x[:__rank])
+    sigma_inv_approx = np.reciprocal(sigma_x)
 
     a_approx = np.linalg.multi_dot(
-        [
-            u_x[:, :__rank].conj().T,
-            x_next,
-            sigma_inv_approx[None] * v_x_t[:__rank, :].conj().T,
-        ]
+        [u_x.conj().T, x_next, sigma_inv_approx[None] * v_x]
     )
 
     return np.linalg.eigvals(a_approx)
@@ -368,11 +257,7 @@ def _varpro_preprocessing(
                   np.ndarray]
     """
 
-    u_r, s_r, v_r_t = np.linalg.svd(data, full_matrices=False)
-    __rank = _compute_rank(s_r, data.shape[0], data.shape[1], rank)
-    u_r = u_r[:, :__rank]
-    s_r = s_r[:__rank]
-    v_r = v_r_t[:__rank].conj().T
+    u_r, s_r, v_r = compute_svd(data, rank)
     data_out = v_r.conj().T * s_r[:, None] if use_proj else data
 
     # trapezoidal derivative approximation

tests/test_varprodmd.py

@@ -3,21 +3,18 @@
 """
 import numpy as np
 import pytest
-from pydmd.utils import compute_rank
+
+from pydmd import VarProDMD
 from pydmd.varprodmd import (
-    _OptimizeHelper,
-    _compute_dmd_rho,
+    OPT_DEF_ARGS,
     _compute_dmd_jac,
-    _compute_rank,
-    _svht,
+    _compute_dmd_rho,
+    _OptimizeHelper,
     compute_varprodmd_any,
     optdmd_predict,
     select_best_samples_fast,
-    OPT_DEF_ARGS,
 )
 
-from pydmd import VarProDMD
-
 
 def signal(x_loc: np.ndarray, time: np.ndarray) -> np.ndarray:
     """
@@ -34,21 +31,6 @@ def signal(x_loc: np.ndarray, time: np.ndarray) -> np.ndarray:
     return __f_1 + __f_2
 
 
-def test_rank():
-    """
-    Test SVHT rank (no duplicate SVD computation is performed).
-    """
-    time = np.linspace(0, 4 * np.pi, 100)
-    x_loc = np.linspace(-10, 10, 1024)
-    __x, __time = np.meshgrid(x_loc, time)
-    z = signal(__x, __time).T
-    s = np.linalg.svd(z, full_matrices=False)[1]
-    assert _compute_rank(s, z.shape[0], z.shape[1], 0) == compute_rank(z, 0)
-    assert _compute_rank(s, z.shape[0], z.shape[1], 4) == 4
-    assert _compute_rank(s, z.shape[0], z.shape[1], 0.8) == compute_rank(z, 0.8)
-    assert _svht(s, z.shape[0], z.shape[1]) == compute_rank(z, 0)
-
-
 def test_varprodmd_rho():
     """
     Unit test for residual vector :math: `\boldsymbol{\rho}`.