Non-orthogonality in SVD solver!

Hi, I am using KrylovKit to compute the SVD of a sparse, asymmetric NXN matrix. The creation of the matrix comes from some previous complicated normalization steps. The SVD is supposed to decompose K as K = U\Lambda V^*, where U, V are both NXL orthogonal matrices.

The orthogonality of the singular vectors is a very essential aspect of SVD. Unfortunately, the SVD computed by KrylovKit does not produce orthogonal singular vectors ! In the code below, I have plotted the inner product matrix A of the first 5 singular vectors stored in matrix Phi.

using Serialization;
using KrylovKit;
using LinearAlgebra;
using RecursiveArrayTools;

file = "Rot_Test_BistochKer_N=10000_Q=0_eps=0.1.dat";
L=201;

K=deserialize(file); print(size(K), "\n");
Lambda, Phi, Gamma, info = KrylovKit.svdsolve( K, L, krylovdim=2*L+1, verbosity=2, tol=1E-5, maxiter=100000 );
L=size(Lambda,1);
if L%2 == 0 L-=1; end
Phi = convert(Array,VectorOfArray(Phi[1:L]));
A = Phi'*Phi; display(A[1:5, 1:5]);

The output is

(10000, 10000)
[ Info: GKL svdsolve in iter 1: 202 values converged, normres = (0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 1.14e-90, 1.21e-140, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 4.92e-314, 0.00e+00, 6.43e-135, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 3.03e-48, 3.56e-238, 0.00e+00, 7.06e-210, 2.09e-197, 0.00e+00, 1.26e-45, 1.22e-191, 0.00e+00, 0.00e+00, 1.44e-119, 0.00e+00, 0.00e+00, 5.57e-240, 4.54e-295, 0.00e+00, 9.61e-50, 0.00e+00, 9.66e-258, 6.45e-309, 0.00e+00, 9.47e-202, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 5.32e-274, 0.00e+00, 0.00e+00, 1.51e-136, 0.00e+00, 0.00e+00, 0.00e+00, 3.74e-15, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 1.72e-15, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 7.90e-193, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 1.01e-213, 5.64e-48, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 7.09e-124, 0.00e+00, 0.00e+00, 7.18e-244, 3.87e-319, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 4.06e-96, 2.27e-46, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 7.08e-133, 6.86e-272, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 3.75e-265, 0.00e+00, 1.48e-203, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 5.33e-201, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 1.49e-244, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 6.61e-46, 0.00e+00, 3.14e-300, 2.02e-310, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 8.16e-283, 2.76e-73, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 2.28e-149, 4.00e-29, 0.00e+00, 4.47e-180, 0.00e+00, 0.00e+00, 0.00e+00, 6.76e-304, 9.51e-226, 0.00e+00, 0.00e+00, 0.00e+00, 8.31e-256, 2.05e-289, 0.00e+00, 8.09e-36, 0.00e+00, 2.78e-107, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 9.70e-199, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 4.42e-231, 0.00e+00, 5.33e-186, 0.00e+00, 0.00e+00)
┌ Info: GKL svdsolve finished after 1 iterations:
│ *  202 singular values converged
│ *  norm of residuals = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1389922138221656e-90, 1.2136856313573638e-140, 0.0, 0.0, 0.0, 0.0, 4.9233869797e-314, 0.0, 6.432194574249315e-135, 0.0, 0.0, 0.0, 0.0, 0.0, 3.033445707684106e-48, 3.563056738382421e-238, 0.0, 7.05736339977992e-210, 2.0888516682420584e-197, 0.0, 1.2581217360698388e-45, 1.220487880475007e-191, 0.0, 0.0, 1.4376102870718524e-119, 0.0, 0.0, 5.567277934850544e-240, 4.540991615018736e-295, 0.0, 9.610688187076416e-50, 0.0, 9.65768779327434e-258, 6.45313901491702e-309, 0.0, 9.472232628461268e-202, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.320549988751628e-274, 0.0, 0.0, 1.5053994469585886e-136, 0.0, 0.0, 0.0, 3.7436841536012125e-15, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.723702476342736e-15, 0.0, 0.0, 0.0, 0.0, 7.903763741365516e-193, 0.0, 0.0, 0.0, 0.0, 1.007372244830932e-213, 5.643244466113458e-48, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.091666692310085e-124, 0.0, 0.0, 7.183284728820492e-244, 3.86903e-319, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.064350229949593e-96, 2.274515101604475e-46, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.075095991722647e-133, 6.861726266719724e-272, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.750383532876986e-265, 0.0, 1.4800348311321002e-203, 0.0, 0.0, 0.0, 0.0, 0.0, 5.331953177945732e-201, 0.0, 0.0, 0.0, 0.0, 0.0, 1.4866236036589676e-244, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.608810497586961e-46, 0.0, 3.1352059855573106e-300, 2.0159614553935e-310, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.1614613388996805e-283, 2.7586069495376143e-73, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.2766045363782646e-149, 4.001460947178077e-29, 0.0, 4.473135387144816e-180, 0.0, 0.0, 0.0, 6.764730684384456e-304, 9.505571784450633e-226, 0.0, 0.0, 0.0, 8.314627790094353e-256, 2.0504203962884657e-289, 0.0, 8.085930143077215e-36, 0.0, 2.778616202307629e-107, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.699566209086634e-199, 0.0, 0.0, 0.0, 0.0, 4.416583944799526e-231, 0.0, 5.325887537657697e-186, 0.0, 0.0, 0.0)
└ *  number of operations = 404
5×5 Array{Float64,2}:
  1.0           0.996786     -1.73472e-18   5.72459e-17  -0.00611994
  0.996786      1.0          -0.00449723    0.00281483    1.51606e-10
 -1.73472e-18  -0.00449723    1.0           5.11743e-17   0.686736
  5.72459e-17   0.00281483    5.11743e-17   1.0          -0.550513
 -0.00611994    1.51606e-10   0.686736     -0.550513      1.0

As you can see, even the first two singular vectors are not orthgonal to each other. The LinearAlgebra has an svd function, but does not handle sparse matrices. Here is a google drive link to the K matrix.

1 Like