Unexpected result from lsqlin

Recently I encountered an unexpecte result from lsqlin. Basically I have a linear system that can be any size, even 1x1. This is actually the case of interest. I attach the code:

s=[7.46355685131195e-06;2.61224489795918e-07;-7.46355685131195e-06;2.61224489795918e-07;0;0;0;0;0;0;0;0;0;0;2.61224489795918e-07;1.21904761904762e-08;-2.61224489795918e-07;6.09523809523810e-09;0;0;0;0;0;0;0;0;0;0;-7.46355685131195e-06;-2.61224489795918e-07;1.49271137026239e-05;0;-7.46355685131195e-06;2.61224489795918e-07;0;0;0;0;0;0;0;0;2.61224489795918e-07;6.09523809523810e-09;0;2.43809523809524e-08;-2.61224489795918e-07;6.09523809523810e-09;0;0;0;0;0;0;0;0;0;0;-7.46355685131195e-06;-2.61224489795918e-07;1.49271137026239e-05;0;-7.46355685131195e-06;2.61224489795918e-07;0;0;0;0;0;0;0;0;2.61224489795918e-07;6.09523809523810e-09;0;2.43809523809524e-08;-2.61224489795918e-07;6.09523809523810e-09;0;0;0;0;0;0;0;0;0;0;-7.46355685131195e-06;-2.61224489795918e-07;1.49271137026239e-05;0;-7.46355685131195e-06;2.61224489795918e-07;0;0;0;0;0;0;0;0;2.61224489795918e-07;6.09523809523810e-09;0;2.43809523809524e-08;-2.61224489795918e-07;6.09523809523810e-09;0;0;0;0;0;0;0;0;0;0;-7.46355685131195e-06;-2.61224489795918e-07;1.87420339044574e-05;8.27690603250169e-08;-1.12784770531454e-05;3.43993550120935e-07;0;0;0;0;0;0;0;0;2.61224489795918e-07;6.09523809523810e-09;8.27690603250169e-08;2.61795472287276e-08;-3.43993550120935e-07;6.99453551912568e-09;0;0;0;0;0;0;0;0;0;0;-1.12784770531454e-05;-3.43993550120935e-07;1.74696244982844e-05;-1.13373307789508e-07;-6.19114744513898e-06;2.30620242331427e-07;0;0;0;0;0;0;0;0;3.43993550120935e-07;6.99453551912568e-09;-1.13373307789508e-07;2.54432097407122e-08;-2.30620242331427e-07;5.72706935123043e-09;0;0;0;0;0;0;0;0;0;0;-6.19114744513898e-06;-2.30620242331427e-07;1.23822948902779e-05;-2.91167575618666e-22;0;0;0;0;0;0;0;0;0;0;2.30620242331427e-07;5.72706935123043e-09;-2.91167575618666e-22;2.29082774049217e-08];
z=[-8287.20136961842;-207.695733121256;-7600.97514253616;-534.391927360151;-5646.81895072033;-842.926473299020;-2749.01799205948;-608.263953041959;-1800.96549474418;-218.233277809774;-1246.61121855062;-222.701572169097;-454.985962888617;-229.489284788134;-207.695733121258;-4.31159651230973;-211.583072054673;-10.5817786869719;-207.469127449374;-19.4084812806167;-205.554052675241;-14.2334094340071;-127.980376432081;-5.82044377019796;-89.3993392954632;-10.0358904762527;-38.6028133285081;-12.8937660172881;-7600.97514253634;-211.583072054673;-6704.83421266254;-532.664018086793;-4247.04164361325;-784.667235831777;-374.396860142614;-575.355545468737;-646.254735224021;-237.767202031333;-565.502081093624;-211.802648689606;-245.788051294455;-215.794095284558;-534.391927360146;-10.5817786869719;-532.664018086794;-24.4250894201275;-533.575869722970;-44.6885886888697;-562.261663238565;-27.6728741709106;-340.294987483632;-3.91648028872448;-233.005638525614;-13.8504678942122;-98.9719393675175;-19.5724723345471;-5646.81895072033;-207.469127449374;-4247.04164361319;-533.575869722965;-874.961989329080;-681.827735041812;4550.00831251650;-524.710100984108;1754.26594854213;-278.456336087023;813.872071787322;-192.126815096214;174.918786479950;-186.341242826303;-842.926473299020;-19.4084812806167;-784.667235831777;-44.6885886888698;-681.827735041812;-73.1296054491854;-560.603026422607;-44.4867666242581;-367.876515600137;-5.36578138470502;-262.187708991285;-13.9626208876454;-116.168335764689;-19.6047309183772;-2749.01799205944;-205.554052675241;-374.396860142644;-562.261663238565;4550.00831251667;-560.603026422608;12023.7350216699;-465.718352496805;5337.97747610125;-338.507025489656;2511.57328179486;-178.101015182985;439.839230323449;-183.023532525701;-608.263953041959;-14.2334094340071;-575.355545468737;-27.6728741709106;-524.710100984101;-44.4867666242582;-465.718352496805;-15.3539182495451;-359.289436100471;17.6126762396493;-285.779177534844;7.60510383356906;-145.366784817360;0.0955479313669716;-1800.96549474418;-127.980376432081;-646.254735224029;-340.294987483632;1754.26594854214;-367.876515600137;5337.97747610143;-359.289436100471;2994.26975751692;-289.642137402146;1958.78945016873;-162.344662349205;811.528922972301;-119.802928873498;-218.233277809774;-5.82044377019797;-237.767202031333;-3.91648028872449;-278.456336087023;-5.36578138470502;-338.507025489655;17.6126762396493;-289.642137402146;35.6258773907798;-248.235934165285;25.5115603160600;-137.758663681854;17.5455698744229;-1246.61121855062;-89.3993392954633;-565.502081093631;-233.005638525614;813.872071787326;-262.187708991284;2511.57328179485;-285.779177534844;1958.78945016873;-248.235934165288;1588.04981195810;-144.928904474322;800.700832420902;-94.0273547749384;-222.701572169097;-10.0358904762527;-211.802648689606;-13.8504678942122;-192.126815096214;-13.9626208876454;-178.101015182985;7.60510383356905;-162.344662349205;25.5115603160600;-144.928904474323;20.8134206727082;-85.8545148794237;15.3928444934107;-454.985962888616;-38.6028133285082;-245.788051294462;-98.9719393675172;174.918786479950;-116.168335764689;439.839230323450;-145.366784817360;811.528922972296;-137.758663681854;800.700832420902;-85.8545148794237;488.887846553582;-53.5142394051840;-229.489284788134;-12.8937660172881;-215.794095284558;-19.5724723345471;-186.341242826303;-19.6047309183772;-183.023532525701;0.0955479313669731;-119.802928873498;17.5455698744229;-94.0273547749366;15.3928444934106;-53.5142394051840;10.5398416193868];
A=s'*s;
b=s'*z;
L=-3e10;
U=1e10;
x1=lsqlin(A,b,[],[],[],[],L,U);
x2=b/A;
disp(A*x1-b)
disp(A*x2-b)

Since the matrix A is actually a scalar, the solution can be computed also as b/A. Accidentally, b/A is also compliant to the lower and upper bound L and U, so I'm wondering why does lsqlin return a different solution. I use matlab 2015b with optimization toolbox 7.3. I thank you very much in advance.

[Torsten]

Did you check whether A and b are the same before and after the call to lsqlin ?

Best wishes

Torsten.

[Roberto]

Yes, they are the same. lsqlin does not change the value of the input. Thanks for your answer.

Roberto

Answer 1

[John D'Errico]

Your problem is that you do not understand what / does, the difference between / and \, or the linear algebra, or what lsqlin does. I'm sorry, but that is true. (I'm not sure how to say that more tactfully, and it would take more time to write. Facts are facts.) Effectively, you have made several mistakes, that all combine into one mess of misunderstandings.

First, s is a column vector, as is z.

The solution that you are using uses something often (sometimes?) called the normal equations, to solve a linear least squares problem. It solves the problem

y = A*x

where A is a known matrix, as is y. (Y is usually a vector of knowns.) The normal equations have you solve this as

x = inv(A'*A)*A'*y

but some people think they are using an improvement is they use a tool like lsqlin here:

x = (A'*A,A'*y);

That is effectively what you have done.

However, a better solution to the problem is simply to use BACKSLASH, thus the \ operator, here:

x = A\y

backslash uses better linear algebra algorithms to solve the problems, NEVER forming the matrix A'*A. The problem is that when you create the matrix A'*A, you square the condition number of A. That makes a bad problem worse, is A is nearly singular.

In this case, the real problem is not that you formed A'*A here, although there are always numerical issues to worry about. The problem is that you think that this is the same as using the FORWARD slash operator!

Given two COLUMN vectors s and z, here is what you did:

A = s'*s
A =
   2.3126e-09

b
b =
     0.069217

When you solved the problem as

x2 = b/A

MATLAB sees that b and A are scalar variables. It simply divides the scalar b by the scalar A, so it is solving the problem as:

/   Slash or right division.
    B/A is the matrix division of A into B, which is roughly the
    same as B*INV(A) , except it is computed in a different way.
    More precisely, B/A = (A'\B')'. See \.

which, in the case that A was a MATRIX, would be a different thing. Luckily, A is a scalar. Slash implicitly solves the problem that would have been posed as

x2*(s'*s) = s'*z

Since x2, (s'*s) and (s'*z) are ALL scalar variables, scalar multiplication commutes, and we can simply use divide here, so forward slash. Lets see what happens though.

format long g
x2 = b/A
x2 =
        29930312.2584127

Here is what BACKSLASH (what you SHOULD have used in the first place) yields

x = s\z
x =
        29930312.2584127

Here is what lsqlin gives, with no constraints:

lsqlin(A,b)
ans =
        29930312.2584127

lsqlin(s,z)
ans =
        29930312.2584127

So they all worked, and gave the same result. LSQLIN has a problem only when you give it bound constraints, due to the fact that a bound constrained problem forces it to use a different algorithm.

L=-3e10;
U=1e10;
x1=lsqlin(A,b,[],[],[],[],L,U)
Optimization terminated: relative function value changing by less
 than sqrt(OPTIONS.FunctionTolerance), and rate of progress is slow.
x1 =
        2.03201411306158

So LSQLIN could have done better, had it simply checked that A and b are scalars, and used a better (read that as simple) algorithm in that case. But LSQLIN was a bit dumb, and used an iterative method, that employs transformations. Those HUGE limits on the solution mess it up terribly.

For example, had I written this:

x1=lsqlin(A,b,[],[],[],[],1e6,1e8)
Optimization terminated: relative function value changing by less
 than OPTIONS.FunctionTolerance.
x1 =
        29930312.2600999

LSQLIN does a decent job, that gave around 9 correct digits or so. It still used an iterative method to solve the problem, when a simple divide, followed by a test to see if the result lay within the bounds would have sufficed.

But your bounds set a huge possible range for the result. So LSQLIN had to transform the problem, then it used an iterative scheme with a tolerance, and got a bit lost. Again, LSQLIN should be more intelligent here. But I think the author of LSQLIN never expected LSQLIN to be used for scalar division, so they never bothered to check for that case.

Let me try lsqlin with more reasonable bounds, to solve a real linear algebra problm as you might have posed it:

x1=lsqlin(s,z,[],[],[],[],1e6,1e8)
Optimization terminated: relative function value changing by less
 than OPTIONS.FunctionTolerance.
x1 =
           29930312.267375

x1=lsqlin(s,z,[],[],[],[],L,U)
Maximum number of iterations exceeded; increase OPTIONS.MaxIterations.
x1 =
                     298.5

These both try to solve the problem

s*x1 = z

where x1 is a scalar unknown variable. The difference is the width of the bounds.

By the way, I'm not sure why you are solving this problem, nor am I sure that s is the vector you REALLY want to be using. Note that roughly half of s is ZERO.

[s,z]
ans =
    7.46355685131195e-06         -8287.20136961842
    2.61224489795918e-07         -207.695733121256
   -7.46355685131195e-06         -7600.97514253616
    2.61224489795918e-07         -534.391927360151
                       0         -5646.81895072033
                       0          -842.92647329902
                       0         -2749.01799205948
                       0         -608.263953041959
                       0         -1800.96549474418
                       0         -218.233277809774
                       0         -1246.61121855062
                       0         -222.701572169097
                       0         -454.985962888617
                       0         -229.489284788134
    2.61224489795918e-07         -207.695733121258
    1.21904761904762e-08         -4.31159651230973
   -2.61224489795918e-07         -211.583072054673
     6.0952380952381e-09         -10.5817786869719
                       0         -207.469127449374
                       0         -19.4084812806167
                       0         -205.554052675241
                       0         -14.2334094340071
                       0         -127.980376432081
                       0         -5.82044377019796
                       0         -89.3993392954632
                       0         -10.0358904762527
                       0         -38.6028133285081
                       0         -12.8937660172881
   -7.46355685131195e-06         -7600.97514253634
   -2.61224489795918e-07         -211.583072054673
    1.49271137026239e-05         -6704.83421266254
                       0         -532.664018086793
   -7.46355685131195e-06         -4247.04164361325
    2.61224489795918e-07         -784.667235831777
                       0         -374.396860142614
                       0         -575.355545468737
...

I'll stop there. So to solve the problem

s*x1 = z

see that roughly half of s is zero, so the residuals for those elements are always in the form of 0*x1 = z(i).

So I would bet that in reality whatever process created the vector s, there are truly some non-zero elements there that you theoretically WANTED to have as non-zero. But hey, what do I know?

Anyway, I expect that you will get upset at my response because I've told you that you don't understand what you are doing. But the facts are that you used the wrong equations to solve a problem, got lucky because scalar multiplication commutes where a matrix multiply would not have done so, then used LSQLIN without understanding what the bounds were doing to you. All of this suggests that you did not really understand what you were doing.

It is true that LSQLIN should be smart enough to see this is a scalar problem, and that they (the authors) could/should have special cased scalar problems without recourse to iterative methods and linear algebra. I'll assume only that nobody ever expected someone would have used lsqlin to divide one number into another number.

[Roberto]

No actually I'm very grateful for your long and detailed answer. Concerning the / vs. \ issue, I know the difference but in the code that I attached I consciously decided to use the / because I knew that the equation was scalar, I'm sorry for not pointing that out. I used lsqlin for scalar division because, as I said in my post, I'd like my code to work for any dimension, also for 1x1 systems, but it's not the typical case, I can have 2x2, 3x3 and so on, which makes the use of lsqlin quite reasonable in my opinion.

Basically my code doesn't work because of the constraints which are too big, I realized that doing some tests after the post. Thank you for your interesting explanation, though. I managed to solve such an issue by "scaling" the system. Basically I solve 10^9*A*x=b and then I multiply the solution by 10^9, so the result is the solution of A*x=b. I don't know if there are numerical drawbacks, apparently in such a way I can set smaller bounds and for what I can tell it works.

Thank you also for your remark about normal equations, I guess I can considerably improve my code without forming A'*A. Nice!

Since you're asking why are there many 0s in s, just reshape s as a 14x14 matrix and you'll see that s has some structure.