_procrustes.py
2.72 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
"""
Solve the orthogonal Procrustes problem.
"""
from __future__ import division, print_function, absolute_import
import numpy as np
from .decomp_svd import svd
__all__ = ['orthogonal_procrustes']
def orthogonal_procrustes(A, B, check_finite=True):
"""
Compute the matrix solution of the orthogonal Procrustes problem.
Given matrices A and B of equal shape, find an orthogonal matrix R
that most closely maps A to B using the algorithm given in [1]_.
Parameters
----------
A : (M, N) array_like
Matrix to be mapped.
B : (M, N) array_like
Target matrix.
check_finite : bool, optional
Whether to check that the input matrices contain only finite numbers.
Disabling may give a performance gain, but may result in problems
(crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns
-------
R : (N, N) ndarray
The matrix solution of the orthogonal Procrustes problem.
Minimizes the Frobenius norm of ``(A @ R) - B``, subject to
``R.T @ R = I``.
scale : float
Sum of the singular values of ``A.T @ B``.
Raises
------
ValueError
If the input array shapes don't match or if check_finite is True and
the arrays contain Inf or NaN.
Notes
-----
Note that unlike higher level Procrustes analyses of spatial data, this
function only uses orthogonal transformations like rotations and
reflections, and it does not use scaling or translation.
.. versionadded:: 0.15.0
References
----------
.. [1] Peter H. Schonemann, "A generalized solution of the orthogonal
Procrustes problem", Psychometrica -- Vol. 31, No. 1, March, 1996.
Examples
--------
>>> from scipy.linalg import orthogonal_procrustes
>>> A = np.array([[ 2, 0, 1], [-2, 0, 0]])
Flip the order of columns and check for the anti-diagonal mapping
>>> R, sca = orthogonal_procrustes(A, np.fliplr(A))
>>> R
array([[-5.34384992e-17, 0.00000000e+00, 1.00000000e+00],
[ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00],
[ 1.00000000e+00, 0.00000000e+00, -7.85941422e-17]])
>>> sca
9.0
"""
if check_finite:
A = np.asarray_chkfinite(A)
B = np.asarray_chkfinite(B)
else:
A = np.asanyarray(A)
B = np.asanyarray(B)
if A.ndim != 2:
raise ValueError('expected ndim to be 2, but observed %s' % A.ndim)
if A.shape != B.shape:
raise ValueError('the shapes of A and B differ (%s vs %s)' % (
A.shape, B.shape))
# Be clever with transposes, with the intention to save memory.
u, w, vt = svd(B.T.dot(A).T)
R = u.dot(vt)
scale = w.sum()
return R, scale