_procrustes.py 2.7 KB

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