[Help-gsl] singular value mismatch between gsl_linalg_SV_decomp() and gsl_linalg_SV_decomp_jacobi()

[Awhan Patnaik]

I have a parameter estimation problem at hand. 'y = d * cons'. Given 'y'
and 'd' I have to estimate 'cons'. In general 'd' is ill-conditioned so
I try to compute pseudo-inverse estimates of 'cons' i.e. 'cons_hat' =
pinv(y)*d. pinv() is Octave's routine to compute the pseudo-inverse. In
real problems I will never know the actual value of 'cons' but in this
test problem I explicitly generate 'cons'.

The pseudo-inverse solution provides a minimum norm estimate so I expect
that the norm of estimated cons i.e. 'cons_hat' should be less than or
equal to 'cons' i.e. norm(cons_hat) <= norm(cons).

d is a 250 by 50 matrix
cons is 50 by 2 matrix
y is a 250 by 2 matrix

I used the gsl_linalg_SV_decomp_jacobi() and gsl_linalg_SV_decomp() to
compute the pseudo-inverse solution but there is a great mismatch
between the two solutions. The singular values of 'd' computed by the
two methods are as follows:

gsl_linalg_SV_decomp()          gsl_linalg_SV_decomp_jacobi()
+3.455928050446208033e+01       +3.455928050446205901e+01
+3.118025106085450560e+01       +3.118025106085450915e+01
+2.876173799945803466e+01       +2.876173799945803822e+01
+2.703614637895373463e+01       +2.703614637895373107e+01
+1.607556212711750021e+01       +1.607556212711749311e+01
+1.028288294325071739e+01       +1.028288294325070851e+01
+8.956977364647210393e+00       +8.956977364647210393e+00
+8.572301064596866027e+00       +8.572301064596869580e+00
+4.416162910061038005e+00       +4.416162910061036229e+00
+4.064036336046190634e+00       +4.064036336046191522e+00
+2.840537542233247059e+00       +2.840537542233246171e+00
+2.590659851516940559e+00       +2.590659851516939227e+00
+2.339470960617429984e+00       +2.339470960617429984e+00
+9.530123960825008789e-01       +9.530123960825004348e-01
+6.801290300942436362e-01       +6.801290300942440803e-01
+5.651059323431912862e-01       +5.651059323431912862e-01
+5.141449512671836253e-01       +5.141449512671836253e-01
+1.647995444010272592e-01       +1.647995444010004196e-01
+1.274166094968336438e-01       +1.274166094968170737e-01
+9.590373999760284929e-02       +9.590373999626004842e-02
+8.827755867423001113e-02       +8.827755867285479174e-02
+2.354167206441644400e-02       +2.354167206115808861e-02
+1.631002297804827125e-02       +1.631002297700456793e-02
+1.283479723818775967e-02       +1.283479723753867471e-02
+1.013834421995763249e-02       +1.013834421883918688e-02
+2.578097221467176743e-03       +2.578097192448009032e-03
+1.585848131462468027e-03       +1.585848121830359331e-03
+1.422983944057921524e-03       +1.422983932711491869e-03
+9.464880833021338492e-04       +9.464879915584590176e-04
+2.159371930184189321e-04       +2.159338662556846097e-04
+1.414251841818201677e-04       +1.414227154760039254e-04
+1.143177118272336419e-04       +1.143166940324424457e-04
+7.092001335167483298e-05       +7.091922842670118216e-05
+1.325147641599603683e-05       +1.299711045877224340e-05
+8.169476585479867420e-06       +6.329950768890471854e-06
+7.059324121096642276e-06       +5.761890661820608645e-06
+3.864854772760345783e-06       +5.503376254333912406e-06
+4.858474059654523714e-07       +3.440395468446005404e-06
+3.526641856963069515e-07       +2.686503336589600847e-06
+2.799341421737882738e-07       +2.218578323076922936e-06
+1.457889938451175779e-07       +2.146523733078234208e-06
+1.905446577885018222e-09       +1.641911328790403411e-06
+2.216073826887362193e-10       +1.641138273866448715e-06
+3.484711650614306555e-12       +1.491486068987434906e-06
+1.893665896724350137e-13       +1.123365345076111080e-06
+8.419908355204420962e-15       +7.714041454208867937e-07
+4.757118843195127556e-16       +7.307018549576065344e-07
+4.467610642465890165e-16       +2.750796959158087445e-07
+3.792345980947965983e-16       +1.669375647399612485e-07
+3.086315108413938371e-16       +8.741960206823821514e-08

Notice how the smallest singular values differ in both cases. One is of
the order of 1e-16 and the other 1e-8. I think the lower 16 singular
values are inconsistent with each other.  GNU Octave seems to agree
with gsl_linalg_SV_decomp() though.

After estimating 'cons_hat' I reconstruct 'y' i.e. y_hat = d * cons_hat

The following is what I observed:

one-norm(cons) = 55.2668 this is the base value,
the minimum norm solution should be smaller or equal to this.

With gsl_linalg_SV_decomp() and singular value cut-off/threshold = 0

one_norm(cons_hat) :: 97.2507 <-- larger than base value
one_norm(y_hat - y) :: 1.67077e-12 <-- good reconstruction

With gsl_linalg_SV_decomp_jacobi() and singular value cut-off/threshold = 0

one_norm(cons_hat_jacobi) :: 4.14156e+06 <-- HUGE ???
one_norm(y_hat - y) :: 2.41439 <-- not good

Next I tried to use a tolerance to get rid of very small singular values.

I first use the MATLAB tolerance MAX(SIZE(A)) * NORM(A) * EPS(class(A)).

With gsl_linalg_SV_decomp() and singular value cut-off/threshold = 1.91843e-12

one_norm(cons_hat) :: 53.5834 <-- smaller than base value, good.
one_norm(y_hat - y) :: 2.20602e-12 <-- good reconstruction


With gsl_linalg_SV_decomp_jacobi() and singular value
cut-off/threshold = 1.91843e-12

one_norm(cons_hat_jacobi) :: 4.14156e+06 <-- still HUGE ???
one_norm(y_hat - y) :: 2.41439 <-- still bad

Next I concentrate on the jacobi version only.

I keep increasing the threshold.


With gsl_linalg_SV_decomp_jacobi() and singular value cut-off/threshold = 1e-6

one_norm(cons_hat_jacobi) :: 14169.1 <-- reduced but not good yet.
one_norm(y_hat - y) :: 0.0507182 <-- improved but not good


With gsl_linalg_SV_decomp_jacobi() and singular value cut-off/threshold = 1e-3

one_norm(cons_hat_jacobi) :: 52.928 <-- now it comes below base value
one_norm(y_hat - y) :: 0.0104453 <-- don't like it.


Well I am plain confused about the situation. The jacobi version is
supposed to be more reliable but in this case how do I decide upon a
suitable threshold. Should I keep on increasing the tolerance value? But
will the solution be a reliable estimate of the parameters?

I have the following questions/comments:


1) Which method should I prefer? gsl_linalg_SV_decomp_jacobi() or
gsl_linalg_SV_decomp()
2) How do I decide the threshold for the jacobi version?
3) If some pinv solution yields a lower norm than the base value
should I simply choose that.
Although in real world problems I will never know the base norm value.
4) Please advise on how to go about solving this problem in some other
way known to you.
5) If you notice any other mistake or omission then please let me know.


Details of the PC I used for the program:
Linux 3.5.4-1-ARCH #1 SMP PREEMPT i686 GNU/Linux
Archlinux 32 bit OS on a 64 bit machine
g++ (GCC) 4.7.1 20120721 (prerelease)
gsl 1.15-2

Compilation command:
g++ main.cpp -lgsl -lgslcblas

I have used a fixed seed of 123UL to generate all the matrices etc.

Below is how the code runs on my PC.

$ ./a.out
d->size1 :: 250
d->size2 :: 50
one_norm(cons) :: 55.2668
cons->size1 :: 50
cons->size2 :: 2

y->size1 :: 250
y->size2 :: 2
one_norm(cons_hat) :: 97.2507
reconstruction_error(cons_hat, d, y) :: 1.67077e-12

one_norm(cons_hat_jacobi) :: 4.14156e+06
reconstruction_error(cons_hat_jacobi, d, y) :: 2.41439

let's use tolerance to discard small singular values

one_norm(cons_hat) :: 53.5834
reconstruction_error(cons_hat, d, y) :: 2.20602e-12

one_norm(cons_hat_jacobi) :: 4.14156e+06
reconstruction_error(cons_hat_jacobi, d, y) :: 2.41439
increase tolerance value for the jacobi method
one_norm(cons_hat_jacobi) :: 14169.1
reconstruction_error(cons_hat_jacobi, d, y) :: 0.0507182

one_norm(cons_hat_jacobi) :: 52.928
reconstruction_error(cons_hat_jacobi, d, y) :: 0.0104453

I have attached the code that generates the data as well as the data
in files themselves.

I will be very grateful for any help received.

Thanks in advance.
Awhan

singular_values.txt
Description: Text document

main.cpp
Description: Text Data

singular_values_jacobi.txt
Description: Text document

y.txt
Description: Text document

d.txt
Description: Text document

cons.txt
Description: Text document