BUG: linalg.cossin: fix for nonsymmetric cases

[Ritvik1sharma]

…omposition correctly. If I run it for a 9x9 matrix with p, q as 3, 6 and 6,3 it causes errors.

Reference issue

What does this implement/fix?

Additional information

[ilayn]

Hi @Ritvik1sharma thank you for taking the time to fix an issue. Could you please elaborate what the issue was and what the solution fixes?

For instance, there is still #19365 open that I just realized I forgot to come back to. Is it possible to replicate the story for this issue too?

If you have a test case that was failing and now does not, we should add it to our test suite.

[melissawm]

Hi @Ritvik1sharma , do you need help to push this forward? Cheers!

[mdhaber]

Actually, I think this just needs tests. The following property-based test attempts to confirm everything the documentation claims about the output. It fails (many times over) in main, but passes with this patch. The patch also fixes gh-19365.

import numpy as np
from scipy import linalg
from scipy.stats import unitary_group

def test_cossin(m, p, q):
    # Generate unitary input
    X = unitary_group.rvs(m)
    np.testing.assert_allclose(X @ X.conj().T, np.eye(m), atol=1e-15)

    # Perform the decomposition
    u, cs, vh = linalg.cossin(X, p=p, q=q, separate=True)
    u1, u2 = u
    v1, v2 = vh
    v1, v2 = v1.conj().T, v2.conj().T

    # "U1, U2, V1, V2 are square orthogonal/unitary matrices
    # of dimensions (p,p), (m-p,m-p), (q,q), and (m-q,m-q) respectively"
    np.testing.assert_allclose(u1 @ u1.conj().T, np.eye(p), atol=1e-13)
    np.testing.assert_allclose(u2 @ u2.conj().T, np.eye(m-p), atol=1e-13)
    np.testing.assert_allclose(v1 @ v1.conj().T, np.eye(q), atol=1e-13)
    np.testing.assert_allclose(v2 @ v2.conj().T, np.eye(m-q), atol=1e-13)

    # "and C and S are (r, r) nonnegative diagonal matrices..."
    C = np.diag(np.cos(cs))
    S = np.diag(np.sin(cs))
    # "...satisfying C^2 + S^2 = I where r = min(p, m-p, q, m-q)."
    r = min(p, m-p, q, m-q)
    np.testing.assert_allclose(C**2 + S**2, np.eye(r))

    # "Moreover, the rank of the identity matrices are
    # min(p, q) - r, min(p, m - q) - r, min(m - p, q) - r,
    # and min(m - p, m - q) - r respectively."
    I11 = np.eye(min(p, q) - r)
    I12 = np.eye(min(p, m - q) - r)
    I21 = np.eye(min(m - p, q) - r)
    I22 = np.eye(min(m - p, m - q) - r)

    # From:
    #                            ┌                   ┐
    #                            │ I  0  0 │ 0  0  0 │
    # ┌           ┐   ┌         ┐│ 0  C  0 │ 0 -S  0 │┌         ┐*
    # │ X11 │ X12 │   │ U1 │    ││ 0  0  0 │ 0  0 -I ││ V1 │    │
    # │ ────┼──── │ = │────┼────││─────────┼─────────││────┼────│
    # │ X21 │ X22 │   │    │ U2 ││ 0  0  0 │ I  0  0 ││    │ V2 │
    # └           ┘   └         ┘│ 0  S  0 │ 0  C  0 │└         ┘
    #                            │ 0  0  I │ 0  0  0 │
    #                            └                   ┘

    # We can see that U and V are block diagonal matrices like so:
    U = linalg.block_diag(u1, u2)
    V = linalg.block_diag(v1, v2)

    # And the center matrix, which we'll call Q here, must be:
    Q11 = np.zeros((u1.shape[1], v1.shape[0]))
    IC11 = linalg.block_diag(I11, C)
    Q11[:IC11.shape[0], :IC11.shape[1]] = IC11

    Q12 = np.zeros((u1.shape[1], v2.shape[0]))
    SI12 = linalg.block_diag(-S, -I12)
    Q12[-SI12.shape[0]:, -SI12.shape[1]:] = SI12

    Q21 = np.zeros((u2.shape[1], v1.shape[0]))
    SI21 = linalg.block_diag(S, I21)
    Q21[-SI21.shape[0]:, -SI21.shape[1]:] = SI21

    Q22 = np.zeros((u2.shape[1], v2.shape[0]))
    IC22 = linalg.block_diag(I22, C)
    Q22[:IC22.shape[0], :IC22.shape[1]] = IC22

    Q = np.block([[Q11, Q12], [Q21, Q22]])

    # Confirm that `cossin` decomposes `X` as shown
    np.testing.assert_allclose(X, U @ Q @ V.conj().T)

    # And check that `separate=False` agrees
    U0, CS0, Vh0 = linalg.cossin(X, p=p, q=q)
    np.testing.assert_allclose(U, U0)
    np.testing.assert_allclose(Q, CS0)
    np.testing.assert_allclose(V, Vh0.conj().T)

for m in range(1, 20):
    for p in range(1, m):
        for q in range(1, m):
            print(p, q)
            test_cossin(m, p, q)

Thanks @Ritvik1sharma!

@ilayn I see that when separate=False, cossin assembles the matrices from the separate pieces in a way that leaves a lot of opportunities for indexing mistakes. Above, I assemble the separate pieces from the block matrices in a way that requires a lot less thinking. Without this patch, the outputs with separate=False don't match mine, and they don't multiply to reproduce the original matrix (in some cases). With this patch, they outputs match mind, and they always multiply to reproduce the original matrix. So I haven't thought carefully about why this works, but it definitely seems to.

Do you want me to add something like the test above to the test suite? There may be some overlap with what is already there, though.

Also, would you like to take a closer look before this merges?

[sfcaracciolo]

Hi everybody,

I want to share with you that I tested the @Ritvik1sharma code too. This fragment https://github.com/sfcaracciolo/cossin_wrapper/blob/master/tests/issue.py works for every q < m now. Before it failed if m-p < q < m.

Thanks!

[ilayn]

I think it is obvious that there is a bug in the indexing but I did not look into it. So any fix is welcome.

[mdhaber]

OK, I'll go ahead and add something like my test above and merge.

[mdhaber]

@sfcaracciolo I think this is about ready, but I might add compute_u and compute_vh into the test for completeness. Otherwise, do the tests look like they are comprehensive and doing what they purport to do (i.e. check all the features about the decomposition claimed in the docs)?

[sfcaracciolo]

@mdhaber I think the tests looks fine, great job!

Just as suggestion, I think it's good to build Q11, Q12, Q21 and Q22 without u1, u2, v1, v2 information (shapes in this case) to keep Q separated of v, u errors. Following this, I would replace,

Q11 = np.zeros((u1.shape[1], v1.shape[0]))
...

by

r = min([p, q, m - p, m - q])
fmin = lambda a, b: min([a, b])-r
fmax = lambda a, b: max([a-b, 0])

Q11 = np.zeros((fmin(p,q) + r + fmax(p,q), fmin(p,q) + r + fmax(q,p)))
...

or something like (adapted of this),

Q11 = sp.sparse.block_diag((
    sp.sparse.identity(fmin(p,q), format='dia', dtype=np.float64),
    C,
    sp.sparse.csc_array((fmax(p,q), fmax(q,p)), dtype=np.float64)
))

[mdhaber]

Just as suggestion, I think it's good to build Q11, Q12, Q21 and Q22 without u1, u2, v1, v2 information (shapes in this case) to keep Q separated of v, u errors.

Thanks! Note that the dimensions of the u and v matrices were already checked. Maybe that protects against the sort of errors you had in mind?

In any case, I was following the properties given in the documentation, and I don't see the shapes of the Q matrices given there explicitly. After this merges, would you be willing to add that information to the documentation and update the test accordingly? I'd be happy to review.

[mdhaber]

@sfcaracciolo I'll go ahead and merge now so this fix makes it into the next release, but please let me know if you'd be willing to add that information to the documentation. I'd be happy to review and merge a PR about that.