Testing for PSDP

procrustes/psdp.py

@@ -23,30 +23,91 @@
 """Positive semidefinite Procrustes Module."""
 
 from math import inf, sqrt
+from typing import Optional
 
 import numpy as np
 from numpy.linalg import multi_dot
-from procrustes.utils import ProcrustesResult
+from procrustes.utils import ProcrustesResult, setup_input_arrays
 import scipy
 from scipy.optimize import minimize
 
 
 __all__ = ["psdp_woodgate"]
 
 
-def psdp_woodgate(a: np.ndarray, b: np.ndarray) -> ProcrustesResult:
+def psdp_woodgate(
+    a: np.ndarray,
+    b: np.ndarray,
+    pad: bool = True,
+    translate: bool = False,
+    scale: bool = False,
+    unpad_col: bool = False,
+    unpad_row: bool = False,
+    check_finite: bool = True,
+    weight: Optional[np.ndarray] = None,
+) -> ProcrustesResult:
     r"""
     Woodgate's algorithm for positive semidefinite Procrustes.
 
+    Given a matrix :math:`\mathbf{A}_{n \times m}` and a reference matrix :math:`\mathbf{B}_{n
+    \times m}`, find the positive semidefinite transformation matrix :math:`\mathbf{P}_{n
+    \times n}` that makes :math:`\mathbf{PA}` as close as possible to :math:`\mathbf{B}`.
+    In other words,
+
+    .. math::
+        \text{PSDP: } min_{\mathbf{P}} \|\mathbf{B}-\mathbf{P}\mathbf{A}\|_{F}^2
+
+    This Procrustes method requires the :math:`\mathbf{A}` and :math:`\mathbf{B}` matrices to
+    have the same shape, which is gauranteed with the default ``pad`` argument for any given
+    :math:`\mathbf{A}` and :math:`\mathbf{B}` matrices. In preparing the :math:`\mathbf{A}` and
+    :math:`\mathbf{B}` matrices, the (optional) order of operations is: **1)** unpad zero
+    rows/columns, **2)** translate the matrices to the origin, **3)** weight entries of
+    :math:`\mathbf{A}`, **4)** scale the matrices to have unit norm, **5)** pad matrices with zero
+    rows/columns so they have the same shape.
+
+    Throughout all methods used for implementing the Woodgate algorithm, the matrices
+    :math:`\mathbf{A}` and :math:`\mathbf{B}` are referred to as :math:`\mathbf{G}` and
+    :math:`\mathbf{F}` respectively, following the nomenclature used in [1]_.
+
     Parameters
     ----------
     a : np.ndarray
-        The matrix to be transformed.
-        This is relabelled to :math:`\mathbf{G}` as in the paper.
+        The matrix :math:`\mathbf{A}` which is to be transformed.
+        This is relabelled to variable g representing the matrix :math:`\mathbf{G}` as
+        in the paper.
 
     b : np.ndarray
-        The target matrix.
-        This is relabellled to :math:`\mathbf{F}` as in the paper.
+        The target matrix :math:`\mathbf{B}`.
+        This is relabelled to variable f representing the matrix :math:`\mathbf{F}` as
+        in the paper.
+
+    pad : bool, optional
+        Add zero rows (at the bottom) and/or columns (to the right-hand side) of matrices
+        :math:`\mathbf{A}` and :math:`\mathbf{B}` so that they have the same shape.
+
+    translate : bool, optional
+        If True, both arrays are centered at origin (columns of the arrays will have mean zero).
+
+    scale : bool, optional
+        If True, both arrays are normalized with respect to the Frobenius norm, i.e.,
+        :math:`\text{Tr}\left[\mathbf{A}^\dagger\mathbf{A}\right] = 1` and
+        :math:`\text{Tr}\left[\mathbf{B}^\dagger\mathbf{B}\right] = 1`.
+
+    unpad_col : bool, optional
+        If True, zero columns (with values less than 1.0e-8) on the right-hand side of the intial
+        :math:`\mathbf{A}` and :math:`\mathbf{B}` matrices are removed.
+
+    unpad_row : bool, optional
+        If True, zero rows (with values less than 1.0e-8) at the bottom of the intial
+        :math:`\mathbf{A}` and :math:`\mathbf{B}` matrices are removed.
+
+    check_finite : bool, optional
+        If True, convert the input to an array, checking for NaNs or Infs.
+
+    weight : ndarray, optional
+        The 1D-array representing the weights of each row of :math:`\mathbf{A}`. This defines the
+        elements of the diagonal matrix :math:`\mathbf{W}` that is multiplied by :math:`\mathbf{A}`
+        matrix, i.e., :math:`\mathbf{A} \rightarrow \mathbf{WA}`.
 
     Returns
     -------
@@ -55,9 +116,9 @@ def psdp_woodgate(a: np.ndarray, b: np.ndarray) -> ProcrustesResult:
 
     Notes
     -----
-    Given :math:`\mathbf{F}, \mathbf{G} \in R^{\ n\times m}`, the woodgate
-    algorithm finds :math:`\mathbf{P} \in S^{\ n\times n}_{\geq}` such
-    that the following is true:
+    Given :math:`\mathbf{F}, \mathbf{G} \in \mathbb{R}^{\ n\times m}`, the woodgate
+    algorithm finds the positive semidefinite matrix :math:`\mathbf{P}_{n \times n}`
+    such that the following is true:
 
     .. math::
         \text{PSDP: } min_{\mathbf{P}} \|\mathbf{F}-\mathbf{P}\mathbf{G}\|
@@ -67,8 +128,8 @@ def psdp_woodgate(a: np.ndarray, b: np.ndarray) -> ProcrustesResult:
     original problem.
 
     .. math::
-        \text{PSDP*: } min_{\mathbf{E} \in R^{\ n\times n}} \|\mathbf{F}
-            - \mathbf{E}^T\mathbf{E}\mathbf{G}\|
+        \text{PSDP*: } min_{\mathbf{E} \in \mathbb{R}^{\ n\times n}}
+            \|\mathbf{F} - \mathbf{E}^T\mathbf{E}\mathbf{G}\|
 
     Now, since all local minimizers of PSDP* are also global minimizers, we
     have :math:`\hat{\mathbf{P}} = \hat{\mathbf{E}}^T\mathbf{E}` where
@@ -94,12 +155,28 @@ def psdp_woodgate(a: np.ndarray, b: np.ndarray) -> ProcrustesResult:
         procrustes". Journal of the American Statistical Association, 93(453),
         584-588.
     """
+    # Check inputs and define the matrices F (matrix to be transformed) and
+    # G (the target matrix).
+    g, f = setup_input_arrays(
+        a,
+        b,
+        unpad_col,
+        unpad_row,
+        pad,
+        translate,
+        scale,
+        check_finite,
+        weight,
+    )
+    if g.shape != f.shape:
+        raise ValueError(
+            f"Shape of A and B does not match: {g.shape} != {f.shape} "
+            "Check pad, unpad_col, and unpad_row arguments."
+        )
 
-    # We define the matrices F, G and Q as in the paper.
+    # We define the matrix Q as in the paper.
     # Our plan is to find matrix P such that, |F - PG| is minimized.
     # Now, |F - PG| = |PG - F| = |PGI - F| = |PAI - F|.
-    f = b
-    g = a
     q = f.dot(g.T) + g.dot(f.T)
 
     # We define the functions L and f as in the paper.
@@ -232,12 +309,12 @@ def _find_gradient(e: np.ndarray, le: np.ndarray, g: np.ndarray) -> np.ndarray:
             = (I\otimes L(\mathbf{E}_i)) v(\mathbf{E}_1)
 
     In the following implementation, the variables d1 and d2 represent
-    :math:`\mathbf{D}_1` and :math:`\mathbf{D}_2`, respectively.
+    :math:`\mathbf{D}_1` and :math:`\mathbf{D}_2`, respectively. Refer to
+    equations (26) and (27) in [1]_ for exact deinitions of the several terms mentioned
+    in this function.
 
     References
     ----------
-    Refer to equations (26) and (27) in [1] for exact deinitions of the
-    several terms mentioned in this function.
     .. [1] Woodgate, K. G. (1996). "A new algorithm for positive
         semidefinite procrustes". Journal of the American Statistical
         Association, 93(453), 584-588.

procrustes/test/test_psdp.py

@@ -0,0 +1,88 @@
+# -*- coding: utf-8 -*-
+# The Procrustes library provides a set of functions for transforming
+# a matrix to make it as similar as possible to a target matrix.
+#
+# Copyright (C) 2017-2022 The QC-Devs Community
+#
+# This file is part of Procrustes.
+#
+# Procrustes is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 3
+# of the License, or (at your option) any later version.
+#
+# Procrustes is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, see <http://www.gnu.org/licenses/>
+#
+# --
+r"""Testings for PSDP (positive semi-definite Procrustes) module."""
+
+import numpy as np
+from numpy.testing import assert_almost_equal
+from procrustes.psdp import psdp_woodgate
+
+
+def test_psdp_identity(n=np.random.randint(50, 100)):
+    r"""Test PSDP with identity matrix."""
+    a = np.eye(n)
+    b = np.eye(n)
+    res = psdp_woodgate(a=a, b=b)
+    s = res["s"]
+    error = res["error"]
+    assert_almost_equal(s, np.eye(n))
+    assert_almost_equal(error, 0.0)
+
+
+def test_psdp_diagonal():
+    r"""Test PSDP with diagonal matrix."""
+    a = np.diag([1, 2, 3, 4])
+    b = np.eye(4)
+    res = psdp_woodgate(a=a, b=b)
+    s = res["s"]
+    error = res["error"]
+    actual_result = np.diag([0.99999, 0.5, 0.33333, 0.25])
+    assert_almost_equal(s, actual_result, decimal=5)
+    assert_almost_equal(error, 0.0, decimal=3)
+
+
+def test_psdp_generic_square():
+    r"""Test PSDP with 2 generic square matrices."""
+    # The example given here is from the original paper.
+    a = np.array([[1, 6, 0], [4, 3, 0], [0, 0, -0.5]])
+    b = np.array([[1, 0, 0], [0, -2, 3], [0, 2, 4]])
+    res = psdp_woodgate(a=a, b=b)
+    s = res["s"]
+    error = res["error"]
+    actual_result = np.array(
+        [
+            [0.22351489, -0.11059539, 0.24342428],
+            [-0.11059539, 0.05472271, -0.12044658],
+            [0.24342428, -0.12044658, 0.26510708],
+        ]
+    )
+    assert_almost_equal(s, actual_result)
+    assert_almost_equal(error, 5.600999068630569, decimal=3)
+
+
+def test_psdp_generic_non_square():
+    r"""Test PSDP with 2 generic non-square matrices."""
+    # The example given here is from the original paper.
+    a = np.array([[5, 1, 6, -1], [3, 2, 0, 2], [2, 4, 3, -3]])
+    b = np.array([[15, 1, 15 - 3, 2 + 5], [10, 5, 6, 3], [-3, 3, -3, -2 + 4]])
+    res = psdp_woodgate(a=a, b=b)
+    s = res["s"]
+    error = res["error"]
+    actual_result = np.array(
+        [
+            [2.57997197, 1.11007896, -1.08770156],
+            [1.11007896, 1.68429863, 0.12829214],
+            [-1.08770156, 0.12829214, 0.75328052],
+        ]
+    )
+    assert_almost_equal(s, actual_result, decimal=5)
+    assert_almost_equal(error, 5.674530443833204, decimal=3)