help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: PetrSt
Subject: Re: Solving A*x=b when A is full rank but numerically rank deficient
Date: Wed, 26 Jun 2013 07:40:57 -0700 (PDT)

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)] &lt; 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. &lt;a
href=&quot;http://gmplib.org/&quot;>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
<http://octave.1599824.n4.nabble.com/file/n4654857/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
<http://octave.1599824.n4.nabble.com/file/n4654857/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
<http://octave.1599824.n4.nabble.com/file/n4654857/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).



--
View this message in context: 
http://octave.1599824.n4.nabble.com/Solving-A-x-b-when-A-is-full-rank-but-numerically-rank-deficient-tp4653604p4654857.html
Sent from the Octave - General mailing list archive at Nabble.com.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]