|
def orthogonal( |
|
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, |
|
lapack_driver: str = "gesvd", |
|
) -> ProcrustesResult: |
|
r"""Perform orthogonal Procrustes. |
|
|
|
Given a matrix :math:`\mathbf{A}_{m \times n}` and a reference matrix :math:`\mathbf{B}_{m |
|
\times n}`, find the orthogonal transformation matrix :math:`\mathbf{Q}_{n |
|
\times n}` that makes :math:`\mathbf{AQ}` as close as possible to :math:`\mathbf{B}`. |
|
In other words, |
|
|
|
.. math:: |
|
\underbrace{\min}_{\left\{\mathbf{Q} | \mathbf{Q}^{-1} = {\mathbf{Q}}^\dagger \right\}} |
|
\|\mathbf{A}\mathbf{Q} - \mathbf{B}\|_{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. |
|
|
|
Parameters |
|
---------- |
|
a : ndarray |
|
The 2D-array :math:`\mathbf{A}` which is going to be transformed. |
|
b : ndarray |
|
The 2D-array :math:`\mathbf{B}` representing the reference matrix. |
|
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}`. |
|
lapack_driver : {'gesvd', 'gesdd'}, optional |
|
Whether to use the more efficient divide-and-conquer approach ('gesdd') or the more robust |
|
general rectangular approach ('gesvd') to compute the singular-value decomposition with |
|
`scipy.linalg.svd`. |
|
|
|
Returns |
|
------- |
|
res : ProcrustesResult |
|
The Procrustes result represented as a class:`utils.ProcrustesResult` object. |
|
|
|
Notes |
|
----- |
|
The optimal orthogonal matrix is obtained by, |
|
|
|
.. math:: |
|
\mathbf{Q}^{\text{opt}} = |
|
\arg \underbrace{\min}_{\left\{\mathbf{Q} \left| {\mathbf{Q}^{-1} = {\mathbf{Q}}^\dagger} |
|
\right. \right\}} \|\mathbf{A}\mathbf{Q} - \mathbf{B}\|_{F}^2 = |
|
\arg \underbrace{\max}_{\left\{\mathbf{Q} \left| {\mathbf{Q}^{-1} = {\mathbf{Q}}^\dagger} |
|
\right. \right\}} \text{Tr}\left[\mathbf{Q^\dagger}\mathbf{A^\dagger}\mathbf{B}\right] |
|
|
|
The solution is obtained using the singular value decomposition (SVD) of the |
|
:math:`\mathbf{A}^\dagger \mathbf{B}` matrix, |
|
|
|
.. math:: |
|
\mathbf{A}^\dagger \mathbf{B} &= \tilde{\mathbf{U}} \tilde{\mathbf{\Sigma}} |
|
\tilde{\mathbf{V}}^{\dagger} \\ |
|
\mathbf{Q}^{\text{opt}} &= \tilde{\mathbf{U}} \tilde{\mathbf{V}}^{\dagger} |
|
|
|
The singular values are always listed in decreasing order, with the smallest singular |
|
value in the bottom-right-hand corner of :math:`\tilde{\mathbf{\Sigma}}`. |
|
|
|
Examples |
|
-------- |
|
>>> import numpy as np |
|
>>> from scipy.stats import ortho_group |
|
>>> from procrustes import orthogonal |
|
>>> a = np.random.rand(5, 3) # random input matrix |
|
>>> q = ortho_group.rvs(3) # random orthogonal transformation |
|
>>> b = np.dot(a, q) + np.random.rand(1, 3) # random target matrix |
|
>>> result = orthogonal(a, b, translate=True, scale=False) |
|
>>> print(result.error) # error (should be zero) |
|
>>> print(result.t) # transformation matrix (same as q) |
|
>>> print(result.new_a) # translated array a |
|
>>> print(result.new_b) # translated array b |
|
|
|
""" |
|
# check inputs |
|
new_a, new_b = setup_input_arrays( |
|
a, |
|
b, |
|
unpad_col, |
|
unpad_row, |
|
pad, |
|
translate, |
|
scale, |
|
check_finite, |
|
weight, |
|
) |
|
if new_a.shape != new_b.shape: |
|
raise ValueError( |
|
f"Shape of A and B does not match: {new_a.shape} != {new_b.shape} " |
|
"Check pad, unpad_col, and unpad_row arguments." |
|
) |
|
# calculate SVD of A.T * B |
|
u, _, vt = scipy.linalg.svd(np.dot(new_a.T, new_b), lapack_driver=lapack_driver) |
|
# compute optimal orthogonal transformation |
|
u_opt = np.dot(u, vt) |
|
# compute one-sided error |
|
error = compute_error(new_a, new_b, u_opt) |
|
|
|
return ProcrustesResult(error=error, new_a=new_a, new_b=new_b, t=u_opt, s=None) |