finding partial solutions from ill-conditioned equations

I am hoping you can point me in the right direction. I have a linear system Ax = b. A is illconditioned. I believe, however that the equation can be used to compute an accurate estimate for some elements of x but probably not all. What algorithm can be used to compute the largest subset of x accurately?

For example, A is 10 by 10 but do to its structure only allows accurate calculation of x_1 and x_2. How to find x_1 and x_2?

Thank you in advance. Greg

Answer 1

[Matt J]

If you know which x(i) are well-determined, you could use lsqminnorm. In this example, we know that only x(1:2) are well-determined, so we can do,

A=rand(10).*[1,1,zeros(1,8)]

A = 10×10

0.7178 0.3339 0 0 0 0 0 0 0 0 0.4073 0.2310 0 0 0 0 0 0 0 0 0.5336 0.7164 0 0 0 0 0 0 0 0 0.4747 0.3006 0 0 0 0 0 0 0 0 0.1531 0.5900 0 0 0 0 0 0 0 0 0.7637 0.7071 0 0 0 0 0 0 0 0 0.4597 0.5534 0 0 0 0 0 0 0 0 0.4173 0.8342 0 0 0 0 0 0 0 0 0.9080 0.9670 0 0 0 0 0 0 0 0 0.0889 0.5736 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

b=3*A(:,1)+2*A(:,2);

x=lsqminnorm(A,b); x=x(1:2)

x = 2×1

3.0000 2.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Matt J]

If you don't know which x(i) belong to the well-determined subset, you would need to do a null space analysis, I think.

A=rand(10); 
A(:,3:end)=randn; %reduce rank
A=A(:,randperm(end));

b=3*A(:,1)+2*A(:,2);

N=null(A)

N = 10×7

0.9016 0.1113 0.2230 -0.0000 -0.0000 0.0000 0.0000 -0.3245 0.8046 0.3498 0.0000 -0.0000 0.0000 0.0000 -0.2443 -0.5613 0.7072 0.0000 0.0000 0.0000 -0.0000 -0.0666 -0.0709 -0.2560 0.8590 -0.2446 -0.0481 -0.0000 -0.0666 -0.0709 -0.2560 -0.2511 -0.2269 0.8279 0.0000 -0.0666 -0.0709 -0.2560 0.0224 0.8941 0.0103 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0666 -0.0709 -0.2560 -0.3151 -0.2113 -0.3950 -0.7071 -0.0666 -0.0709 -0.2560 -0.3151 -0.2113 -0.3950 0.7071 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Here, we see that the rows 7 and 10 of N are numerically 0. This tells us that the well-determined components are x(7) and x(10), because no matter what null vector you add to x, x(7) and x(10) will not change. So, you can now do,

x=lsqminnorm(A,b);
x=x([7,10]);

Answer 2

[John D'Errico]

The largest subset? It may be impossible to compute ANY subset accurately. You may be able to only compute a linear combination of the variables with any accuracy. For example:

A = rand(10,1)*[1 1] + randn(10,2)*1e-16;
b = randn(10,1);

If you look at A, you will see the first and second columns are virtually identical, down to about the 16th decimal place.

format long g
A

A = 10×2

0.737647957843521 0.737647957843521 0.650311923800328 0.650311923800327 0.540775422404753 0.540775422404753 0.362989155847819 0.362989155847819 0.00692348027389468 0.0069234802738949 0.63415258209194 0.63415258209194 0.446695664232614 0.446695664232614 0.824746413661878 0.824746413661878 0.0963008279107997 0.0963008279107997 0.572767760139755 0.572767760139755

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now, suppose you want to estimate the vector X=[x1;x2], such that

A*X == b

using linear least squares? Essentially, you cannot estimate either of the variables x1 or x2 well. We can try this:

Xslash = A\b

Warning: Rank deficient, rank = 1, tol = 3.852654e-15.

Xslash = 2×1

0.243288144512921 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Xsvd = lsqminnorm(A,b)

Xsvd = 2×1

0.12164407225646 0.12164407225646

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

There simply is no information content in your data to estimate either x1 or x2 well, while ignoring the other. The only information content in this data lies effectively in the direction x1+x2. And you can do that reasonably well.

Unfortunately, this will often be the case. We can learn what information is lacking by use of svd, but null will do too.

null(A)

ans = 2×1

-0.707106781186547 0.707106781186548

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

And that tells me what you cannot learn is any information about the DIFFERENCE of the two unknowns, as that is where the nullspace lies.

It would have been as easy to generate different data that would fail in some different direction. The problem is, much of the time when you have an ill-conditioned system, there will be no single variables you can estimate truly well, and the others will be less well known. This is easy enough to build such a classic example.

x = rand(100,1);
A = x.^(0:20); % A simple high order polynomial model, constant term first
y = randn(100,1);

A\y

Warning: Rank deficient, rank = 20, tol = 2.220446e-13.

ans = 21×1

1.0e+00 * 0.258903042248013 -50.3895984299905 7467.76406077624 -382744.525828343 9910253.85219477 -154116845.978847 1578723812.50204 -11268772090.3706 58122676388.7107 -221744603150.771

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now, we can try using a tool like lsqminnorm, and I suppose you can say the lowest order terms (constant term, linear term) are at least somewhat better estimated than the rest. The others get worse and worse as you go down the line, but there is no simple dividing line where I would decide which of the unknowns are well-estimated.

lsqminnorm(A,y)

ans = 21×1

1.0e+00 * 0.195007696572898 -29.4286998687226 5539.90134218198 -299865.970766404 7877135.25975838 -122379890.240141 1240835329.82936 -8696578254.88423 43649594747.4703 -160124664408.22

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[John D'Errico]

@James McDaniel

I think what you misunderstand is it is almost never the case where one variable is well defined, and others are not. Almost always it is that some LINEAR combination of the unknowns is poorly defined.

For example, it may be that there is no information available about the linear combination a-2*b. What this means is you can choose ANY value of a, and this will imply the value of b as a/2. In that case, while pinv or lsqminnorm will give you AN answer, it will not tell you the slightest bit of help on what any specific variable must be.