How do I keep the order of vectors after symmetric orthogonalization with Scipy?

Question

scipy.linalg.orth orthogonalizes matrices with SVD. After orthogonalization, the vectors are put in the descending order of singular values. However, for some reasons, I need the orthogonalized vectors to be as similar to the original ones as possible, both in shapes and in orders.

For example, for the matrix,

>>> import scipy.linalg as sl
>>> import numpy as np
>>> A=np.array([[5,0,0],[0,2,0],[0,0,4]])
>>> A
array([[5, 0, 0],
       [0, 2, 0],
       [0, 0, 4]])

instead of getting the orthogonalized result with the descending order of the singular values 5, 4, 2,

>>> sl.orth(A)
array([[1., 0., 0.],
       [0., 0., 1.],
       [0., 1., 0.]])

I want the result to be

1 0 0
0 1 0
0 0 1

which keeps the order of the original vectors. Is there any way I can realize this? Thanks for any help!

Answer 1

You can't directly ask SciPy to do this. SciPy internally uses a LAPACK function called gesdd to do this, which already does this sorting, and does not have a way to disable it.

However, what you could do instead is to compare the cosine similarity of every orthogonal vector to every input vector, and sort the orthogonal vectors by similarity to the input vectors.

Example:

import scipy.linalg
from scipy.spatial.distance import cdist
import numpy as np


def orth_preserving_order(A):
    orthogonals_unordered = scipy.linalg.orth(A)
    # Compute distance matrix
    distance = cdist(A, orthogonals_unordered, metric='cosine')
    # Convert cosine distance into cosine similarity
    similarity_matrix = 1 - distance
    # Consider vectors which are pointing in opposite direction to be "similar"
    similarity_matrix = np.abs(similarity_matrix)
    orthogonal_order_idx = []
    for i in range(A.shape[0]):
        vec_to_match = A[i]
        # Pick element in orthogonals_unordered which is most similar to this row of A
        best_vec_idx = similarity_matrix[i].argmax()
        orthogonal_order_idx.append(best_vec_idx)
        # Prevent this vector from being picked again
        similarity_matrix[:, best_vec_idx] = -np.inf

    # Sanity check - did we put any orthogonal vector in twice?
    assert len(set(orthogonal_order_idx)) == len(orthogonal_order_idx), "duplicate orthogonal_order_idx"
    # Did we use every vector?
    assert len(orthogonal_order_idx) == A.shape[0], "missing vector in orthogonal_order_idx"
    # Reorder orthogonals
    orthogonal_reordered = orthogonals_unordered[orthogonal_order_idx]
    return orthogonal_reordered


A=np.array([[5,0,0],[0,2,0],[0,0,4]])
print(orth_preserving_order(A))

Output:

[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]