Quick Start

Orthogonal Procrustes

Given matrix \(\mathbf{A}_{m \times n}\) and a reference matrix \(\mathbf{B}_{m\times n}\), find the orthogonal transformation matrix \(\mathbf{Q}_{n\times n}\) that makes \(\mathbf{AQ}\) as close as possible to \(\mathbf{B}\), i.e.,

\begin{equation} \underbrace{\min}_{\left\{\mathbf{Q} | \mathbf{Q}^{-1} = {\mathbf{Q}}^\dagger \right\}} \|\mathbf{A}\mathbf{Q} - \mathbf{B}\|_{F}^2 \end{equation}

The code block below gives an example of the orthogonal Procrustes problem for random matrices \(\mathbf{A}\) and \(\mathbf{B}\). Here, matrix \(\mathbf{B}\) is constructed by shifting an orthogonal transformation of matrix \(\mathbf{A}\), so the matrices can be perfectly matched. As is the case with all Procrustes flavours, the user can specify whether the matrices should be translated (so that both are centered at origin) and/or scaled (so that both are normalized to unity with respect to the Frobenius norm). In addition, the other optional arguments (not appearing in the code-block below) specify whether the zero columns (on the right-hand side) and rows (at the bottom) should be removed prior to transformation.

[1]:

import numpy as np
from scipy.stats import ortho_group
from procrustes import orthogonal

# random input 10x7 matrix A
a = np.random.rand(10, 7)

# random orthogonal 7x7 matrix T (acting as Q)
t = ortho_group.rvs(7)

# target matrix B (which is a shifted AT)
b = np.dot(a, t) + np.random.rand(1, 7)

# orthogonal Procrustes analysis with translation
result = orthogonal(a, b, scale=True, translate=True)

# compute transformed matrix A (i.e., A x Q)
aq = np.dot(result.new_a, result.t)

# display Procrustes results
print("Procrustes Error = ", result.error)
print("\nDoes the obtained transformation match variable t? ", np.allclose(t, result.t))
print("Does AQ and B matrices match?", np.allclose(aq, result.new_b))

print("Transformation Matrix T = ")
print(result.t)
print("")
print("Matrix A (after translation and scaling) = ")
print(result.new_a)
print("")
print("Matrix AQ = ")
print(aq)
print("")
print("Matrix B (after translation and scaling) = ")
print(result.new_b)

Procrustes Error =  2.4289236009459004e-30

Does the obtained transformation match variable t?  True
Does AQ and B matrices match? True
Transformation Matrix T =
[[-0.05485823 -0.21987681  0.71054591 -0.25910115  0.54815253 -0.11057903
  -0.25285756]
 [-0.52504088  0.57359009  0.08125893 -0.38628239 -0.16691215  0.36377965
  -0.28162755]
 [-0.13427477 -0.35050672  0.5114134  -0.07667139 -0.71683955  0.04718091
   0.27496942]
 [ 0.19650197  0.65967521  0.42739622  0.33295066  0.08143714 -0.10148903
   0.46449961]
 [ 0.45237142 -0.06057057 -0.04420777 -0.48171449  0.15466827  0.61863151
   0.38866554]
 [-0.0206257   0.12956287 -0.16552502 -0.64649468 -0.03687109 -0.66740601
   0.30107121]
 [ 0.67794881  0.21015927  0.12230059 -0.13005245 -0.35481889 -0.12155183
  -0.56892541]]

Matrix A (after translation and scaling) =
[[-0.00055723 -0.01377505 -0.08991782 -0.16477186  0.14750515  0.00704046
  -0.20861997]
 [ 0.13709756 -0.02863268 -0.04751636  0.07709446 -0.14385904  0.10799371
  -0.0032753 ]
 [-0.16853924 -0.17772675  0.06648402 -0.16374724  0.04464136  0.09435058
   0.03478612]
 [ 0.11811112 -0.01507476 -0.06673255  0.139217   -0.0976698  -0.11659321
  -0.03462358]
 [-0.15903118 -0.01484409 -0.16695146  0.11677978 -0.19205803  0.01023694
   0.10550636]
 [ 0.02382673 -0.01263136  0.05182648 -0.10599329 -0.08805356 -0.07338563
   0.00990156]
 [-0.16720248 -0.04508176  0.03784015  0.10676988  0.04590812  0.00066062
  -0.03174471]
 [-0.06525965  0.14948796  0.07651695  0.02690549 -0.13166814 -0.24761947
  -0.19722096]
 [ 0.12908648  0.16851992 -0.01048313 -0.21501857  0.21002468  0.11819247
   0.14823474]
 [ 0.15246788 -0.01024143  0.14893372  0.18276434  0.20522923  0.09912353
   0.17705575]]

Matrix AQ =
[[-0.08789303 -0.13682353 -0.15112392 -0.09097679  0.14960896  0.1194413
   0.08089827]
 [-0.04048379  0.04296137  0.09182028  0.00475755  0.09519934 -0.19631547
  -0.02539345]
 [ 0.10328743 -0.17937652 -0.1835174  -0.03432131 -0.13263092 -0.06584305
   0.06085581]
 [-0.02749881  0.06414488  0.12745368  0.15361799  0.12791064 -0.01422005
  -0.03266849]
 [ 0.04631795  0.19713971 -0.12997567  0.17079912 -0.02302644 -0.14601308
  -0.07885918]
 [-0.05406852 -0.10266454  0.01435798  0.04801396 -0.04504059 -0.00072564
  -0.09940151]
 [ 0.04797432  0.05870918 -0.06350457  0.07497127 -0.08421826  0.02485656
   0.1510763 ]
 [-0.26805587  0.02546819  0.03909608  0.21141601 -0.05464003  0.17025779
  -0.00558307]
 [ 0.05666167 -0.03614484 -0.00256305 -0.36619058  0.00816506  0.1013846
  -0.14997974]
 [ 0.22375865  0.06658609  0.25795658 -0.17208723 -0.04132777  0.00717706
   0.09905505]]

Matrix B (after translation and scaling) =
[[-0.08789303 -0.13682353 -0.15112392 -0.09097679  0.14960896  0.1194413
   0.08089827]
 [-0.04048379  0.04296137  0.09182028  0.00475755  0.09519934 -0.19631547
  -0.02539345]
 [ 0.10328743 -0.17937652 -0.1835174  -0.03432131 -0.13263092 -0.06584305
   0.06085581]
 [-0.02749881  0.06414488  0.12745368  0.15361799  0.12791064 -0.01422005
  -0.03266849]
 [ 0.04631795  0.19713971 -0.12997567  0.17079912 -0.02302644 -0.14601308
  -0.07885918]
 [-0.05406852 -0.10266454  0.01435798  0.04801396 -0.04504059 -0.00072564
  -0.09940151]
 [ 0.04797432  0.05870918 -0.06350457  0.07497127 -0.08421826  0.02485656
   0.1510763 ]
 [-0.26805587  0.02546819  0.03909608  0.21141601 -0.05464003  0.17025779
  -0.00558307]
 [ 0.05666167 -0.03614484 -0.00256305 -0.36619058  0.00816506  0.1013846
  -0.14997974]
 [ 0.22375865  0.06658609  0.25795658 -0.17208723 -0.04132777  0.00717706
   0.09905505]]

Accessing Scaling- and Translation-Related Information

In order to access the scaling factor of the initial two input matrices, we can use _scale_array function in utils module.

[2]:

# accessing the scaling factor when preparing two input matrices
from procrustes.utils import _scale_array

# array_a is scaled to match array_b's norm; after scaling the norm of $a$ equals the norm of $b$
scaled_a, scaling_factor = _scale_array(a, b)
# print(f"The scaled matrix $a$ is {scaled_a}.")
print(f"The scaling factor is {scaling_factor}.")

The scaling factor is 0.8623410931853571.

The user may also want to get the translation matrix or the translated matrix. This can be achieved by using _translate_array function in the utils module.

[3]:

# accessing the translation matrix
from procrustes.utils import _translate_array

# array_a is translated to centroid of array_b
# after translation, the centroid of $a$ will match the centroid of $b$
trans_array_a, centroid = _translate_array(a, b)
# print(f"The translated array $a$ is {trans_array_a}.")
print(f"The centroid is {centroid}.")

The centroid is [-0.36907112 -1.20438156 -0.25281312 -0.3612893  -1.00236417 -0.150441
 -0.93179808].

The corresponding file can be obtained from: