Solving A*x=b when A is full rank but numerically rank deficient

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Solving A*x=b when A is full rank but numerically rank deficient

CdeMills
Hello,

I have to solve a problem implying number of molecules in various states. The issue is the the smallest number is around 10^4 and the biggest 10^34. And the interesting phenomena occur with the smallest fraction.

Here is a pseudo-problem generation:
# generate a symmetric non-singular matrix
A=randn(5,5); A=A.'*A;
 AA=A*kron(logspace(34, 4, 5), ones(5,1));
X=rand(5,1); B=AA*X;
[AA\B X]

The solution clearly originates from a minimum-norm solver: only one component of X is OK. The main issue is that A singular values span the original range of 10^4 to 10^34. Any idea on how to solve this kind of problem overcoming the numerical resolution issue ?

Regards

Pascal

Reply | Threaded
Open this post in threaded view
|

Re: Solving A*x=b when A is full rank but numerically rank deficient

PetrSt
Banned User
If your aim is to get x == X, where x denotes numerical solution of AA*x=B for AA and B defined in your example, then it is impossible.

Firstly some notation, which I'll use:
Ax=b; A=[aij]; x=[xj]; b=[bi];
s( ) = absolute error
r( ) = relative error
A\b = arbitrary method solvin linear system
~ means "similar" or "comparable to", e.g. x~y is meant as x is of the similar order as y
AA, B, X = variables defined in your example

Obtained x=AA\B (x~=X) is correct solution at the used computational precision.  It can be shown by
   x = AA\B; sB = AA*x - B,
where sB is typically in order of 1e18 (note that result of (1e34+5e17)==1e34 is true).
In case of double the precision should be 16 digits, so for the example round(log10(A(i,:)) == [34, 27, 19, 12, 4] and round(log10(bi)) = 34,  Ax = b+s(b) gives s(b)~1e18 and any aij*xj<s(b) becomes insignificant.
It restricts also the number of valid digits for all aij and xj. If error level (variance of Ax and/or b) of Ax=b is given as s(b)^2 = sum_j{[aij*s(xj)]^2+[s(aij)*xj]^2} and numbers are represented by a decadic decomposition as y=sum_k(y(k)*10^n(k)), its obvious that any [aij(k)*10^n(k)].*[xj(q)*10^n(q)] < 1e18 turns to be negligible in sense of equality Ax==b.
For orders of magnitude of aij [1e34,1e27,1e19,1e12,1e4] the highest obtainable precision on x is [1e-16,1e-9,1e-1,1e6,1e12]. It's clear that x(4) and x(5) could be never estimated correctly, if expected value is in between 0 and 1.
For last precise digit of ai1 with order of 1e34/1e16=1e18 the highest precision applicable to x1 is 1e0, thus even for the coefficient of the highest order not all valid digits of the coefficient are significant. For requested precision of x1 say 1e-4 the last significant digit of ai1 is at place of about 1e22.
The above conclusions about obtainable precision of x can be illustrated by a single step of Jacobi's point-relaxation iterative method with X as the initial guess.
   D = diag(diag(AA));
   xe = inv(D)*((D-AA)*X + B);
The difference between X and xe is in perfect line with the maximal expectable precisions of [1e-16,1e-9,1e-1,1e6,1e12].

The first idea how to reach the aim x==X is to use higher precision of calculation. It can be found some tools to do that (e.g. <a href="http://gmplib.org/">http://gmplib.org/), but since I'm not familiar with those and it was particularly what attracted my attantion I wrote something myself (see uploaded files below). Calculation with variable precision will be noted infP(). Surprisingly x=infP(AA\B) is close to AA\B instead of X (note infP(A\b) uses simple Gauss's elimination). The reason is in the fact that B generated by AA*X is already loaded with the error of order 1e18 as can be seen by comparison of B=AA*X and infP(B=AA*X). Solving the system x=infP(AA\infP(B=AA*X)) gives the desired result x==X (in double prec.).

The conclusion can be made that the solution x=AA\B is correct for given computational precision, so that A*x = b + s(b) with the maximal precision (in sence of orders) of s(xj) = log10(aij)/s(bi). To be sure that all estimated xj are grater than their s(xj), we can quess the needed calculation precision by demanding the lowest acceptable precision of x to be grater or equal to the round-off error of calculation, i.e. s(xn)>=max_j(s(aij*xj))~s(bi), where n denotes the index of the component with coefficient of the lowest order (in example n=5 and s(x5)>=s(ai1*x1)~s(bi)). It can be expressed as
   (x1*s(a1))^2 + (s(x1)*a1)^2 <= (xn*s(an))^2 + (s(xn)*an)^2
   (x1*ar*r())^2 + (x1*a1*r())^2 <= (xn*an*r())^2 + (s(xn)*an)^2
         r() = r(aj) = r(xj) = r(b) = relative precision of calculation, e.g. for double r()=1e-16
         it could be assumed that x1*ar*r() == xn*an*r(), then
         s(xn) is left unchanged, since it represents the lowest requested precision and not a round-off error
                 given by computation
   x1*a1*r() <= s(xn)*an
   x1*r() = s(x1) <= s(xn) * an/a1
   if we assume x1~1, then r() <= s(xn) * an/a1
So if we want x5 to be precise up to 1e-4, we heve to use r() <= 1e-4 * 1e4/1e34 = 1e-34. It is only rough guess, therefore smaller r() is desirable. The important point is that with this precision it must be estimated/calculated/measured/known also the righthand side (b) as shown above. I feel this condition as the major obstacle in practical problems. If it is not your case, then I hope it helps you.

File with code for variable prec. calculations:
    operation.m
    It isn't "all input errors proof" and also it isn't optimized anyhow, so please be patient during evaluation.
    I wrote it because of curiosity and not for a real use. It's advisable to use some existing package or tool for variable precision calculations in case of solving your real problems.

File with demo / checking of operation.m performance / correctnes:
   OperPerformance.m
   Note the last example of estimation of e, how the conversion to double changes a precise number.

File dealing with your example:
   AAx_B_demo.m
   The order of printouts follows text of this reply. Actualy this is the only script you need to run (assuming the operation.m is saved on the path).
Reply | Threaded
Open this post in threaded view
|

Re: Solving A*x=b when A is full rank but numerically rank deficient

CdeMills
Sorry to contradict you. I went one step further:

A=randn(5,5); A=A.'*A;
AA=A*kron(logspace(34, 4, 5), ones(5,1)); AAmp = AA+mp(0);
X=rand(5,1); B=AA*X;Bmp=AAmp*X;
AAinv=inv(AAmp.'*AAmp);

double([AA\B AAinv*(AAmp.'*Bmp) X])
ans =

   9.4146e-02   9.4146e-02   9.4146e-02
   2.9772e-09   7.1297e-01   7.1297e-01
   9.4146e-17   7.0114e-01   7.0114e-01
   2.9772e-24   1.5793e-01   1.5793e-01
   9.4146e-32   9.4513e-01   9.4513e-01

As you may notice, ALL the solutions are correctly computed in the second case. 'mp' is a class I'm busy working on, implementing multi-precision arithmetic. I have to explicitelly invert the design matrix, as the left division operator still has some issues.

Regards

Pascal
Reply | Threaded
Open this post in threaded view
|

Re: Solving A*x=b when A is full rank but numerically rank deficient

withaar
In reply to this post by CdeMills
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Solving A*x=b when A is full rank but numerically rank deficient

PetrSt
Banned User
In reply to this post by CdeMills
CdeMills wrote
Sorry to contradict you. I went one step further:
I don't see any contradiction to what I wrote. In cotrary I think it proves it.
As I understand to new lines:
AAmp = AA+mp(0);   % convert AA to mp-class
Bmp=AAmp*X;         % generate B as mp-class by precise calculation of AA*X
AAinv=inv(AAmp.'*AAmp);
x =  AAinv*(AAmp.'*Bmp);  % solve the linear system in least squares sence as x=inv(A'*A)*(A'*B), where all variables are of mp-class, i.e. variable (infinite) precision calculation.

CdeMills wrote
As you may notice, ALL the solutions are correctly computed in the second case.
It agrees with the last example ("infP calculation with infP B") in the "Ax_B_demo.m" send before, where even the Gauss's elimination gives desirable result x==X for square A.
What I wanted to point out is necessity of knowing A and B precisely up to more than 35 digits to achieve all x with minimal precision around 1e-4. Bmp=AAmp*X surely satisfies this condition. If you are dealing with some theoretical modeling etc. and such precision is accesible for your design matrix and righthand side vector, then you don't need to read any further, since the variable prec. calculation solves the problem.
I'm dealing with problems where A is represented by experimental conditions or literature values and B is some measured state variable, all obtained with some experimental error. In this case the idea of precisions higher than let say 15 orders is from the dreamworld. If it threats to your real calculations, you should examine the behaviour of solution with
Bmp = B + mp(0);    % conversion to mp-class of B obtined by double prec. computation of AA*X.
Bmp = AAmp*X + randn(size(X)).*abserr;   % abserr = absolute error expectable for B (abserr~eps().*B~1e18 for the line above).

Regards
Petr
Reply | Threaded
Open this post in threaded view
|

Re: Solving A*x=b when A is full rank but numerically rank deficient

CdeMills
PetrSt wrote
CdeMills wrote
Sorry to contradict you. I went one step further:
I don't see any contradiction to what I wrote. In cotrary I think it proves it.
As I understand to new lines:
AAmp = AA+mp(0);   % convert AA to mp-class
Bmp=AAmp*X;         % generate B as mp-class by precise calculation of AA*X
AAinv=inv(AAmp.'*AAmp);
x =  AAinv*(AAmp.'*Bmp);  % solve the linear system in least squares sence as x=inv(A'*A)*(A'*B), where all variables are of mp-class, i.e. variable (infinite) precision calculation.
In this case, mp is a class with a precision of 250 bits instead of the usual 64 bits.
PetrSt wrote
CdeMills wrote
As you may notice, ALL the solutions are correctly computed in the second case.
It agrees with the last example ("infP calculation with infP B") in the "Ax_B_demo.m" send before, where even the Gauss's elimination gives desirable result x==X for square A.
What I wanted to point out is necessity of knowing A and B precisely up to more than 35 digits to achieve all x with minimal precision around 1e-4. Bmp=AAmp*X surely satisfies this condition. If you are dealing with some theoretical modeling etc. and such precision is accesible for your design matrix and righthand side vector, then you don't need to read any further, since the variable prec. calculation solves the problem.
I'm dealing with problems where A is represented by experimental conditions or literature values and B is some measured state variable, all obtained with some experimental error. In this case the idea of precisions higher than let say 15 orders is from the dreamworld. If it threats to your real calculations, you should examine the behaviour of solution with
Bmp = B + mp(0);    % conversion to mp-class of B obtined by double prec. computation of AA*X.
Bmp = AAmp*X + randn(size(X)).*abserr;   % abserr = absolute error expectable for B (abserr~eps().*B~1e18 for the line above).

Regards
Petr
I do not pretend I measure with THAT level of precision. I just have to deal with equations mixing quantities whose RELATIVE accuracy are the same (measurement at 1% say) but whose orders of magnitude  cover a range of 10^30. With only 64 bits of resolution, most of the terms in the intermediate computations simply vanish due to truncation error. Stated in another way, the numerical accuracy times the biggest value is greater than the smaller values. With 250 bits of resolution, round-off error accumulates at a level smaller than the absolute uncertainty on the smallest values, so the relative accuracy of the final result is at least as good as the relative accuracy of the input data.

Regards

Pascal