Re: Leasqr matrix singular to machine precision

[Archambault Fabien]

Archambault Fabien a écrit :
> Hi,
>
> I am working on octave since a few days because I need to fit a 
> potential on Lennard-Jones parameters (for those who understand).
>
> So I created octave files that can be run in bash.
>
> I think my equations are well written but if possible someone can check.
>
> The script is :
>
> #!/usr/bin/octave -qf
> source "octave.script"
> x = [ 1.367; 0.224500; 1.768200; -0.120; -0.046000; -0.152100 ];
> r(1,1) = RTEMP(1,14);
> r(1,2) = RTEMP(2,14);
> r(1,3) = RTEMP(3,14);
> y(1) = DELTA(14);
> r(2,1) = RTEMP(1,1);
> r(2,2) = RTEMP(2,1);
> r(2,3) = RTEMP(3,1);
> y(2) = DELTA(1);
> r(3,1) = RTEMP(1,5);
> r(3,2) = RTEMP(2,5);
> r(3,3) = RTEMP(3,5);
> y(3) = DELTA(5);
> r(4,1) = RTEMP(1,13);
> r(4,2) = RTEMP(2,13);
> r(4,3) = RTEMP(3,13);
> y(4) = DELTA(13);
> r(5,1) = RTEMP(1,9);
> r(5,2) = RTEMP(2,9);
> r(5,3) = RTEMP(3,9);
> y(5) = DELTA(9);
> r(6,1) = RTEMP(1,12);
> r(6,2) = RTEMP(2,12);
> r(6,3) = RTEMP(3,12);
> y(6) = DELTA(12);
> function y = f(r,x)
>     y = sqrt(x(4)*x(5)) * ( ((x(1)/2+x(2)/2) ./ r(:,1)).^12 - 2 * 
> ((x(1)/2+x(2)/2) ./ r(:,1)).^6 ) \
>           + sqrt(x(4)*x(5)) * ( ((x(1)/2+x(2)/2) ./ r(:,2)).^12 - 2 * 
> ((x(1)/2+x(2)/2) ./ r(:,2)).^6 ) \
>           + sqrt(x(4)*x(6)) * ( ((x(1)/2+x(3)/2) ./ r(:,3)).^12 - 2 * 
> ((x(1)/2+x(3)/2) ./ r(:,3)).^6 );
> endfunction
> [f1,result,kvg,iter] = leasqr(r,y,x,"f")
>
> See in attach file for octave.script.
>
> When I launch the script the message is :
>
> warning: in 
> /usr/share/octave/2.1.73/site/m/octave-forge/optim/leasqr.m near line 
> 299, column 5:
> warning: division by zero
> warning: division by zero
> warning: inverse: matrix singular to machine precision, rcond = 
> 6.93366e-19
> f1 =
>
>   -4.8789e-04
>   -7.2189e-02
>   -2.5896e-02
>   -9.8740e-04
>   -1.0229e-02
>   -2.1931e-03
>
> result =
>
>   1.367000
>   0.224500
>   1.768200
>  -0.120000
>  -0.046000
>  -0.152100
>
> kvg = 1
> iter = 1
>
> The version of octave is 2.1.73 on Mandriva 2007.1. The script leasqr 
> is taken from octave-forge (in rpm files and tested with script 
> available online).
>
> I have searched also on the mailing list but did not find any clue for 
> this problem.
>
> Thanks in advance for your answer,
> Fabien Archambault
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www.cae.wisc.edu/mailman/listinfo/help-octave
Hi again,

I have worked on my scripts since this morning and found some mistakes 
on my distances so it was impossible to fit with the function.

So I modified it :
#!/usr/bin/octave -qf
source "octave.script"
x = [ 1.367; 0.224500; 1.768200; -0.120; -0.046000; -0.152100 ];
r(1,1) = RTEMP(1,8);
r(1,2) = RTEMP(2,8);
r(1,3) = RTEMP(3,8);
y(1) = DELTA(8);
r(2,1) = RTEMP(1,7);
r(2,2) = RTEMP(2,7);
r(2,3) = RTEMP(3,7);
y(2) = DELTA(7);
r(3,1) = RTEMP(1,13);
r(3,2) = RTEMP(2,13);
r(3,3) = RTEMP(3,13);
y(3) = DELTA(13);
r(4,1) = RTEMP(1,2);
r(4,2) = RTEMP(2,2);
r(4,3) = RTEMP(3,2);
y(4) = DELTA(2);
r(5,1) = RTEMP(1,12);
r(5,2) = RTEMP(2,12);
r(5,3) = RTEMP(3,12);
y(5) = DELTA(12);
r(6,1) = RTEMP(1,5);
r(6,2) = RTEMP(2,5);
r(6,3) = RTEMP(3,5);
y(6) = DELTA(5);
function y = f(r,x)
       y = \
       + sqrt( x(4) * x(5) )*( ( (x(1)+x(2)) ./ r(:,1) ).^12 - 2 * ( 
(x(1)+x(2)) ./ r(:,1) ).^6) \
       + sqrt( x(4) * x(5) )*( ( (x(1)+x(2)) ./ r(:,2) ).^12 - 2 * ( 
(x(1)+x(2)) ./ r(:,2) ).^6) \
       + sqrt( x(4) * x(6) )*( ( (x(1)+x(3)) ./ r(:,3) ).^12 - 2 * ( 
(x(1)+x(3)) ./ r(:,3) ).^6) \
;
endfunction
[f1,result,kvg,iter] = leasqr(r,y,x,"f")

And in attach file the new distances.

I have the same error :
warning: in /usr/share/octave/2.1.73/site/m/octave-forge/optim/leasqr.m 
near line 299, column 5:
warning: division by zero
warning: division by zero
warning: inverse: matrix singular to machine precision, rcond = 6.1851e-23
f1 =

   1.66389 +  1.61324i
   3.05701 +  1.87630i
  -0.06294 +  0.22328i
  60.99300 +  2.95108i
  -0.10449 +  0.42837i
   9.93354 +  2.48499i

result =

    1.43557 +   0.00269i
   -0.36752 -   1.95093i
    1.84242 -   0.00268i
   -0.13855 +   0.00670i
  204.68952 +   9.96074i
   -0.23220 +   0.00222i

kvg = 1
iter = 9

Now it makes iterations but the solutions provided are imaginary !

Do I have to add options like the derivate of the function.

Thanks again for any answer.

--
Fabien Archambault
Equipe de dynamique des assemblages membranaires
Unité Mixte de Recherches CNRS UHP 7565
Université Henri-Poincaré, Nancy I BP 239,
54506 Vandoeuvre-lès-Nancy, cedex France

Tél : 03.83.68.43.96

GROUPANOM1 = [ 20 ];
GROUPAX1 = [ 0 ];
GROUPAY1 = [ 0 ];
GROUPAZ1 = [ 0 ];
GROUPANOM2 = [ 20 ];
GROUPAX2 = [ 0 ];
GROUPAY2 = [ 0 ];
GROUPAZ2 = [ 0 ];
GROUPANOM3 = [ 20 ];
GROUPAX3 = [ 0 ];
GROUPAY3 = [ 0 ];
GROUPAZ3 = [ 0 ];
GROUPANOM4 = [ 20 ];
GROUPAX4 = [ 0 ];
GROUPAY4 = [ 0 ];
GROUPAZ4 = [ 0 ];
GROUPANOM5 = [ 20 ];
GROUPAX5 = [ 0 ];
GROUPAY5 = [ 0 ];
GROUPAZ5 = [ 0 ];
GROUPANOM6 = [ 20 ];
GROUPAX6 = [ 0 ];
GROUPAY6 = [ 0 ];
GROUPAZ6 = [ 0 ];
GROUPANOM7 = [ 20 ];
GROUPAX7 = [ 0 ];
GROUPAY7 = [ 0 ];
GROUPAZ7 = [ 0 ];
GROUPANOM8 = [ 20 ];
GROUPAX8 = [ 0 ];
GROUPAY8 = [ 0 ];
GROUPAZ8 = [ 0 ];
GROUPANOM9 = [ 20 ];
GROUPAX9 = [ 0 ];
GROUPAY9 = [ 0 ];
GROUPAZ9 = [ 0 ];
GROUPANOM10 = [ 20 ];
GROUPAX10 = [ 0 ];
GROUPAY10 = [ 0 ];
GROUPAZ10 = [ 0 ];
GROUPANOM11 = [ 20 ];
GROUPAX11 = [ 0 ];
GROUPAY11 = [ 0 ];
GROUPAZ11 = [ 0 ];
GROUPANOM12 = [ 20 ];
GROUPAX12 = [ 0 ];
GROUPAY12 = [ 0 ];
GROUPAZ12 = [ 0 ];
GROUPANOM13 = [ 20 ];
GROUPAX13 = [ 0 ];
GROUPAY13 = [ 0 ];
GROUPAZ13 = [ 0 ];
GROUPANOM14 = [ 20 ];
GROUPAX14 = [ 0 ];
GROUPAY14 = [ 0 ];
GROUPAZ14 = [ 0 ];
GROUPANOM15 = [ 20 ];
GROUPAX15 = [ 0 ];
GROUPAY15 = [ 0 ];
GROUPAZ15 = [ 0 ];
GROUPBNOM1 = [ 1,1,8 ];
GROUPBX1 = [ -.756358,.756358,0 ];
GROUPBY1 = [ 0,0,0 ];
GROUPBZ1 = [ -2.588115,-2.588115,-1.900000 ];
GROUPBNOM2 = [ 1,1,8 ];
GROUPBX2 = [ -.756358,.756358,0 ];
GROUPBY2 = [ 0,0,0 ];
GROUPBZ2 = [ 
-2.68811500000000000000,-2.68811500000000000000,-2.00000000000000000000 ];
GROUPBNOM3 = [ 1,1,8 ];
GROUPBX3 = [ -.756358,.756358,0 ];
GROUPBY3 = [ 0,0,0 ];
GROUPBZ3 = [ 
-2.78811500000000000000,-2.78811500000000000000,-2.10000000000000000000 ];
GROUPBNOM4 = [ 1,1,8 ];
GROUPBX4 = [ -.756358,.756358,0 ];
GROUPBY4 = [ 0,0,0 ];
GROUPBZ4 = [ 
-2.88811500000000000000,-2.88811500000000000000,-2.20000000000000000000 ];
GROUPBNOM5 = [ 1,1,8 ];
GROUPBX5 = [ -.756358,.756358,0 ];
GROUPBY5 = [ 0,0,0 ];
GROUPBZ5 = [ 
-2.98811500000000000000,-2.98811500000000000000,-2.30000000000000000000 ];
GROUPBNOM6 = [ 1,1,8 ];
GROUPBX6 = [ -.756358,.756358,0 ];
GROUPBY6 = [ 0,0,0 ];
GROUPBZ6 = [ 
-3.08811500000000000000,-3.08811500000000000000,-2.40000000000000000000 ];
GROUPBNOM7 = [ 1,1,8 ];
GROUPBX7 = [ -.756358,.756358,0 ];
GROUPBY7 = [ 0,0,0 ];
GROUPBZ7 = [ 
-3.18811500000000000000,-3.18811500000000000000,-2.50000000000000000000 ];
GROUPBNOM8 = [ 1,1,8 ];
GROUPBX8 = [ -.756358,.756358,0 ];
GROUPBY8 = [ 0,0,0 ];
GROUPBZ8 = [ 
-3.28811500000000000000,-3.28811500000000000000,-2.60000000000000000000 ];
GROUPBNOM9 = [ 1,1,8 ];
GROUPBX9 = [ -.756358,.756358,0 ];
GROUPBY9 = [ 0,0,0 ];
GROUPBZ9 = [ 
-3.38811500000000000000,-3.38811500000000000000,-2.70000000000000000000 ];
GROUPBNOM10 = [ 1,1,8 ];
GROUPBX10 = [ -.756358,.756358,0 ];
GROUPBY10 = [ 0,0,0 ];
GROUPBZ10 = [ 
-3.48811500000000000000,-3.48811500000000000000,-2.80000000000000000000 ];
GROUPBNOM11 = [ 1,1,8 ];
GROUPBX11 = [ -.756358,.756358,0 ];
GROUPBY11 = [ 0,0,0 ];
GROUPBZ11 = [ 
-3.58811500000000000000,-3.58811500000000000000,-2.90000000000000000000 ];
GROUPBNOM12 = [ 1,1,8 ];
GROUPBX12 = [ -.756358,.756358,0 ];
GROUPBY12 = [ 0,0,0 ];
GROUPBZ12 = [ 
-4.18811500000000000000,-4.18811500000000000000,-3.50000000000000000000 ];
GROUPBNOM13 = [ 1,1,8 ];
GROUPBX13 = [ -.756358,.756358,0 ];
GROUPBY13 = [ 0,0,0 ];
GROUPBZ13 = [ 
-4.68811500000000000000,-4.68811500000000000000,-4.00000000000000000000 ];
GROUPBNOM14 = [ 1,1,8 ];
GROUPBX14 = [ -.756358,.756358,0 ];
GROUPBY14 = [ 0,0,0 ];
GROUPBZ14 = [ 
-5.18811500000000000000,-5.18811500000000000000,-4.50000000000000000000 ];
GROUPBNOM15 = [ 1,1,8 ];
GROUPBX15 = [ -.756358,.756358,0 ];
GROUPBY15 = [ 0,0,0 ];
GROUPBZ15 = [ 
-5.68811500000000000000,-5.68811500000000000000,-5.00000000000000000000 ];
DELTA = [ 
100.03,65.9,43.63,28.9,19.07,12.48,8.05,5.07,3.08,1.76,0.9,-0.34,-0.2,-0.09,-0.04
 ];
RTEMP = 
[sqrt((GROUPAX1(1)-1*GROUPBX1(1))^2+(GROUPAY1(1)-1*GROUPBY1(1))^2+(GROUPAZ1(1)-1*GROUPBZ1(1))^2),sqrt((GROUPAX2(1)-1*GROUPBX2(1))^2+(GROUPAY2(1)-1*GROUPBY2(1))^2+(GROUPAZ2(1)-1*GROUPBZ2(1))^2),sqrt((GROUPAX3(1)-1*GROUPBX3(1))^2+(GROUPAY3(1)-1*GROUPBY3(1))^2+(GROUPAZ3(1)-1*GROUPBZ3(1))^2),sqrt((GROUPAX4(1)-1*GROUPBX4(1))^2+(GROUPAY4(1)-1*GROUPBY4(1))^2+(GROUPAZ4(1)-1*GROUPBZ4(1))^2),sqrt((GROUPAX5(1)-1*GROUPBX5(1))^2+(GROUPAY5(1)-1*GROUPBY5(1))^2+(GROUPAZ5(1)-1*GROUPBZ5(1))^2),sqrt((GROUPAX6(1)-1*GROUPBX6(1))^2+(GROUPAY6(1)-1*GROUPBY6(1))^2+(GROUPAZ6(1)-1*GROUPBZ6(1))^2),sqrt((GROUPAX7(1)-1*GROUPBX7(1))^2+(GROUPAY7(1)-1*GROUPBY7(1))^2+(GROUPAZ7(1)-1*GROUPBZ7(1))^2),sqrt((GROUPAX8(1)-1*GROUPBX8(1))^2+(GROUPAY8(1)-1*GROUPBY8(1))^2+(GROUPAZ8(1)-1*GROUPBZ8(1))^2),sqrt((GROUPAX9(1)-1*GROUPBX9(1))^2+(GROUPAY9(1)-1*GROUPBY9(1))^2+(GROUPAZ9(1)-1*GROUPBZ9(1))^2),sqrt((GROUPAX10(1)-1*GROUPBX10(1))^2+(GROUPAY10(1)-1*GROUPBY10(1))^2+(GROUPAZ10(1)-1*GROUPBZ10(1))^2),sqrt((GROUPAX11(1)-1*GROUPBX11(1))^2+(GROUPAY11(1)-1*GROUPBY11(1))^2+(GROUPAZ11(1)-1*GROUPBZ11(1))^2),sqrt((GROUPAX12(1)-1*GROUPBX12(1))^2+(GROUPAY12(1)-1*GROUPBY12(1))^2+(GROUPAZ12(1)-1*GROUPBZ12(1))^2),sqrt((GROUPAX13(1)-1*GROUPBX13(1))^2+(GROUPAY13(1)-1*GROUPBY13(1))^2+(GROUPAZ13(1)-1*GROUPBZ13(1))^2),sqrt((GROUPAX14(1)-1*GROUPBX14(1))^2+(GROUPAY14(1)-1*GROUPBY14(1))^2+(GROUPAZ14(1)-1*GROUPBZ14(1))^2),sqrt((GROUPAX15(1)-1*GROUPBX15(1))^2+(GROUPAY15(1)-1*GROUPBY15(1))^2+(GROUPAZ15(1)-1*GROUPBZ15(1))^2);sqrt((GROUPAX1(1)-1*GROUPBX1(2))^2+(GROUPAY1(1)-1*GROUPBY1(2))^2+(GROUPAZ1(1)-1*GROUPBZ1(2))^2),sqrt((GROUPAX2(1)-1*GROUPBX2(2))^2+(GROUPAY2(1)-1*GROUPBY2(2))^2+(GROUPAZ2(1)-1*GROUPBZ2(2))^2),sqrt((GROUPAX3(1)-1*GROUPBX3(2))^2+(GROUPAY3(1)-1*GROUPBY3(2))^2+(GROUPAZ3(1)-1*GROUPBZ3(2))^2),sqrt((GROUPAX4(1)-1*GROUPBX4(2))^2+(GROUPAY4(1)-1*GROUPBY4(2))^2+(GROUPAZ4(1)-1*GROUPBZ4(2))^2),sqrt((GROUPAX5(1)-1*GROUPBX5(2))^2+(GROUPAY5(1)-1*GROUPBY5(2))^2+(GROUPAZ5(1)-1*GROUPBZ5(2))^2),sqrt((GROUPAX6(1)-1*GROUPBX6(2))^2+(GROUPAY6(1)-1*GROUPBY6(2))^2+(GROUPAZ6(1)-1*GROUPBZ6(2))^2),sqrt((GROUPAX7(1)-1*GROUPBX7(2))^2+(GROUPAY7(1)-1*GROUPBY7(2))^2+(GROUPAZ7(1)-1*GROUPBZ7(2))^2),sqrt((GROUPAX8(1)-1*GROUPBX8(2))^2+(GROUPAY8(1)-1*GROUPBY8(2))^2+(GROUPAZ8(1)-1*GROUPBZ8(2))^2),sqrt((GROUPAX9(1)-1*GROUPBX9(2))^2+(GROUPAY9(1)-1*GROUPBY9(2))^2+(GROUPAZ9(1)-1*GROUPBZ9(2))^2),sqrt((GROUPAX10(1)-1*GROUPBX10(2))^2+(GROUPAY10(1)-1*GROUPBY10(2))^2+(GROUPAZ10(1)-1*GROUPBZ10(2))^2),sqrt((GROUPAX11(1)-1*GROUPBX11(2))^2+(GROUPAY11(1)-1*GROUPBY11(2))^2+(GROUPAZ11(1)-1*GROUPBZ11(2))^2),sqrt((GROUPAX12(1)-1*GROUPBX12(2))^2+(GROUPAY12(1)-1*GROUPBY12(2))^2+(GROUPAZ12(1)-1*GROUPBZ12(2))^2),sqrt((GROUPAX13(1)-1*GROUPBX13(2))^2+(GROUPAY13(1)-1*GROUPBY13(2))^2+(GROUPAZ13(1)-1*GROUPBZ13(2))^2),sqrt((GROUPAX14(1)-1*GROUPBX14(2))^2+(GROUPAY14(1)-1*GROUPBY14(2))^2+(GROUPAZ14(1)-1*GROUPBZ14(2))^2),sqrt((GROUPAX15(1)-1*GROUPBX15(2))^2+(GROUPAY15(1)-1*GROUPBY15(2))^2+(GROUPAZ15(1)-1*GROUPBZ15(2))^2);sqrt((GROUPAX1(1)-1*GROUPBX1(3))^2+(GROUPAY1(1)-1*GROUPBY1(3))^2+(GROUPAZ1(1)-1*GROUPBZ1(3))^2),sqrt((GROUPAX2(1)-1*GROUPBX2(3))^2+(GROUPAY2(1)-1*GROUPBY2(3))^2+(GROUPAZ2(1)-1*GROUPBZ2(3))^2),sqrt((GROUPAX3(1)-1*GROUPBX3(3))^2+(GROUPAY3(1)-1*GROUPBY3(3))^2+(GROUPAZ3(1)-1*GROUPBZ3(3))^2),sqrt((GROUPAX4(1)-1*GROUPBX4(3))^2+(GROUPAY4(1)-1*GROUPBY4(3))^2+(GROUPAZ4(1)-1*GROUPBZ4(3))^2),sqrt((GROUPAX5(1)-1*GROUPBX5(3))^2+(GROUPAY5(1)-1*GROUPBY5(3))^2+(GROUPAZ5(1)-1*GROUPBZ5(3))^2),sqrt((GROUPAX6(1)-1*GROUPBX6(3))^2+(GROUPAY6(1)-1*GROUPBY6(3))^2+(GROUPAZ6(1)-1*GROUPBZ6(3))^2),sqrt((GROUPAX7(1)-1*GROUPBX7(3))^2+(GROUPAY7(1)-1*GROUPBY7(3))^2+(GROUPAZ7(1)-1*GROUPBZ7(3))^2),sqrt((GROUPAX8(1)-1*GROUPBX8(3))^2+(GROUPAY8(1)-1*GROUPBY8(3))^2+(GROUPAZ8(1)-1*GROUPBZ8(3))^2),sqrt((GROUPAX9(1)-1*GROUPBX9(3))^2+(GROUPAY9(1)-1*GROUPBY9(3))^2+(GROUPAZ9(1)-1*GROUPBZ9(3))^2),sqrt((GROUPAX10(1)-1*GROUPBX10(3))^2+(GROUPAY10(1)-1*GROUPBY10(3))^2+(GROUPAZ10(1)-1*GROUPBZ10(3))^2),sqrt((GROUPAX11(1)-1*GROUPBX11(3))^2+(GROUPAY11(1)-1*GROUPBY11(3))^2+(GROUPAZ11(1)-1*GROUPBZ11(3))^2),sqrt((GROUPAX12(1)-1*GROUPBX12(3))^2+(GROUPAY12(1)-1*GROUPBY12(3))^2+(GROUPAZ12(1)-1*GROUPBZ12(3))^2),sqrt((GROUPAX13(1)-1*GROUPBX13(3))^2+(GROUPAY13(1)-1*GROUPBY13(3))^2+(GROUPAZ13(1)-1*GROUPBZ13(3))^2),sqrt((GROUPAX14(1)-1*GROUPBX14(3))^2+(GROUPAY14(1)-1*GROUPBY14(3))^2+(GROUPAZ14(1)-1*GROUPBZ14(3))^2),sqrt((GROUPAX15(1)-1*GROUPBX15(3))^2+(GROUPAY15(1)-1*GROUPBY15(3))^2+(GROUPAZ15(1)-1*GROUPBZ15(3))^2)
 ];