SVD not working as expected with complex matrix in scipy?-CodePudding

I'm trying to find approximate nonzero solutions of M@x = 0 using SVD in scipy where M is a complex-valued 4x4 matrix.

First a toy example:

M = np.array([
    [1,1,1,1],
    [1,0,1,1],
    [1,-1,0,0],
    [0,0,0,1e-10]
])

U, s, Vh = scipy.linalg.svd(M)
print(s) # [2.57554368e 00 1.49380718e 00 3.67579714e-01 7.07106781e-11]
print(Vh[-1]) # [ 0.00000000e 00  2.77555756e-16 -7.07106781e-01  7.07106781e-01]
print(np.linalg.norm( M@Vh[-1] )) # 7.07106781193738e-11

So in this case, the smallest (last) value in s is very small, and corresponding last column Vh[-1] is the approximate solution to M@x=0, where M@Vh[-1] is also very small, roughly same order as s[-1].

Now the real example which doesn't work the same way:

M = np.array([[ 1.68572560e-01-3.98053448e-02j,  5.61165939e-01-1.22638499e-01j,
         3.39625823e-02-1.16216469e 00j,  2.65140034e-06-4.10296457e-06j],
       [ 4.17991622e-01 1.33504182e-02j, -4.79190633e-01-2.08562169e-01j,
         4.87429517e-01 3.68070222e-01j, -3.63710538e-05 6.43912577e-06j],
       [-2.18353842e 06-4.20344071e 05j, -2.52806647e 06-2.08794519e 05j,
        -2.01808847e 06-1.96246695e 06j, -5.77147300e-01-3.12598394e 00j],
       [-3.03044160e 05-6.45842521e 04j, -6.85879183e 05 2.07045473e 05j,
         6.14194217e 04-1.28864668e 04j, -7.08794838e 00 9.70230041e 00j]])
U, s, Vh = scipy.linalg.svd(M)
print(s) # [4.42615634e 06 5.70600901e 05 4.68468171e-01 5.21600592e-13]
print(Vh[-1]) # [-5.35883825e-05 0.00000000e 00j  3.74712739e-05-9.89288566e-06j  4.03111556e-06 7.59306578e-06j -8.20834667e-01 5.71165865e-01j]
print(np.linalg.norm( M@Vh[-1] )) # 35.950705194666476

What's going on here? s[-1] is very small, so M@x should have a solution in principle, but Vh[-1] doesn't look like a solution. Is this an issue with M and Vh being complex numbers? A numerical stability/accuracy issue? Something else?

I'd really like to figure out what x would give M@x with roughly the same order of magnitude as s[-1], please let me know any way to solve this.

CodePudding user response：

You forgot the conjugate transpose

The decomposition given by SVD is np.allclose(M, U @ np.diag(s) @ Vh), if s[-1] is small it means that the last column of U @ np.diag(s) ~ M @ np.inv(Vh) ~ M @ Vh.T.conj(). So you can find the use

M @ Vh[-1].T.conj() # [-7.77136331e-14-3.74441041e-13j,
                    #  4.67810503e-14 3.45797987e-13j,
                    #  -2.84217094e-14-1.06581410e-14j,
                    #  7.10542736e-15 3.10862447e-15j]