Last active
August 29, 2015 14:20
-
-
Save bheklilr/032d8278d1b8fef9d97d to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| def matrix_sqrt(A): | |
| '''Solve for B in A = B*B''' | |
| # B can be decomposed into V * S * V.I where S is diagonal and V is a column | |
| # matrix of the eigenvectors (provided it meets certain criteria, but that | |
| # probably won't happen). | |
| V = np.matrix(np.linalg.eig(A)[1]) # Returns eigenvalues and eigenvectors | |
| VI = V.I | |
| # Then | |
| # (V * S * V.I) * (V * S * V.I) = V * S * S * V.I = A | |
| # S * S = V.I * A * V | |
| # And S * S can be solved by taking the square root of each element, since | |
| # S is diagonal. Thus | |
| SS = VI * A * V | |
| # Explicitly converting to diagonal vector, then converting back. Helps with | |
| # floating point rounding errors. If it makes you uncomfortable then uncomment | |
| # assert np.allclose(SS, np.diagflat(np.diag(SS))) | |
| S = np.matrix(np.diagflat(np.sqrt(np.diag(SS)))) | |
| # Now we can calculate B as V * S * V.I | |
| B = V * S * VI | |
| # And in case you don't believe me, here's another assert that all this worked | |
| # assert np.allclose(A, B * B) | |
| return B | |
| A = np.matrix(np.random.random((4, 4)) + 1j * np.random.random((4, 4))) | |
| print(matrix_sqrt(A)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment