LinearSolve precision quirks (Was: Re: Lyapunov Equation)
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg878] LinearSolve precision quirks (Was: Re: Lyapunov Equation)
- From: danl (Daniel Lichtblau)
- Date: Thu, 27 Apr 1995 01:50:41 -0400
- Organization: Wolfram Research, Inc.
In article <3n9t8r$poj at news0.cybernetics.net> Allan Hayes
<hay at haystack.demon.co.uk> writes:
>
> I have been looking over the recent postings about the solution of
> the Liapunov equation
> p . a + Transpose[ a ] . p == -q (1)
> for p; where a is a matrix, p is a symmetric matrix and q is
> symmetric positive definite matrix.
>
> In mg[729 ] Paul Rubin uses Solve to find the entries in p (he has
> q on the right instead of -q);
> In mg[713] Bernd Cebulski obtains a nice formula for the solution.
>
> I wondered whether it would be better to substitute LinearSolve for
> Solve in Paul's solution.
> The investigation turned up some questions that I had not expected
> to have to deal with
> (a) an apparent bug in LinearSolve;
> (b) a query about reliability of solutions.;
> (c) can we allow the precision estimates use in arbitrary
> precision arithmetic to be influenced by the algorithm being used.
>
> I apologise for the length of the following account
>
> The following picks out the working part of Pauls's code ( with q
> changed to -q)
>
> In[1]:=
> LpSolve[a_,q_]:=
> Flatten[
> p/.
> Solve[
> p . a + Transpose[a] . p == -q,
> Union[Flatten[p]]
> ],
> 1
> ]
>
> Here is my LinearSolve code
>
> In[2]:=
> LpSolve2[a_,q_] :=
> p/.Thread[(Union@@p) ->
> LinearSolve[
> Outer[Coefficient,
> Flatten[p.a + Transpose[a] . p],
> Union@@p
> ],
> Flatten[-q]
> ] ]
>
> Test inputs are generated by
>
> In[3]:=
> test[n_] :=
> ( SeedRandom[1];
> p = Array[ Sort[e[##]]&,{n,n}];
> q = IdentityMatrix[n]; (*could be arbitrary positive definite*)
> a = Table[ 2 Random[] - .5, {i, n}, {j, n} ]
> )
> (the technique for constructing the symmetric matrix p might be of
> independent use)
>
> Try the codes with n=2.
>
> In[4]:= test[2];
>
> (a) BUG IN LINEAR SOLVE?
> In[5]:=
> LpSolve2[a,q];
> LinearSolve::nosol: Linear equation encountered
> which has no solution.
> This is presumably due to some test used in LinearSolve.
> Work round: use arbitrary precision arithmetic.
> In[6]:=
> LpSolve2@@SetPrecision[{a,q}, $MachinePrecision +1]
> Out[6]=
> {{-0.1485960799538594, -0.3531050582315298},
> {-0.3531050582315298, 0.357157867314548}}
>
> The following table shows that
> * LpSolve2 is faster than LpSolve for n = 2 to 9
> * the advantage decreases,
> * the precision of the solution, steadily decreases.(for n
> = 9 it is zero).
> The last is to be expected from arbitrary precision arithmetic,
> which gives only the precision that can be conservatively
> guaranteed.
>
> In[7]:=
> data =
> TableForm[
> Table[
> test[n];
> {n, Timing[s[n] = LpSolve[a,q]][[1]],
> Timing[s2[n] =
> LpSolve2@@SetPrecision[{a,q},
> $MachinePrecision +1]
> ][[1]],
> Precision[s2[n]]
> },
> {n,2,9}
> ],
> TableHeadings ->
> { {},
> {"n", "LpSolve\ntime\n", "LpSolve2\ntime", "Precision
> of \nLpSolve2"}
> }
> ]
>
> Out[7]//TableForm=
> n LpSolve LpSolve2 Precision of
> time time LpSolve2
>
> 2 0.116667 Second 0.0666667 Second 15
>
> 3 0.45 Second 0.183333 Second 14
>
> 4 1.3 Second 0.633333 Second 13
>
> 5 2.96667 Second 1.45 Second 10
>
> 6 6.01667 Second 3.3 Second 8
>
> 7 10.7667 Second 7.18333 Second 7
>
> 8 18.4667 Second 12.9333 Second 3
>
> 9 30.4167 Second 26.0833 Second 0
>
> The same degree of loss of precision shows up if we use arbitrary
> precision with LpSolve ,also the speed goes down quite a bit..
>
The LpSolve2 code works fine in our development version. Note, however,
that the mix of overdetermined systems and machine precision will often
fail; the only reason it works at all is that the redundancy in the eqns
is in the form of identical eqns (rather than nontrivial linear
combinations that reduce to 0==0). Small bugs in version 2.2 probably
account for the fact that these obviously consistant eqns failed.
> (b) RELIABILITY OF SOLUTIONS?
> If this loss of precision represents the true precision of the
> algorithm then prsumably the machine arithmetic we used in the table
> above for LpSolve gives a n unreliable answer.
> However, when we substitute s[9] for p on the right side of (1) we get
>
> <omitted >
> (c) IS THE REDUCTION OF PRECISION WHEN USING ARBITRARY PRECISION
> ARITHMETIC EXCESSIVE IN THESE CALCULATIONS, AND, IF SO, WHAT CAN BE
> DONE ABOUT IT?
>
>
> <omitted >
The estimated precision loss is often excessive. If you can call
LinearSolve with a square system of full rank, you are likely to do much
better, at least in our development version. This is because we will now
make better use of special case code. For example, if I change LpSolve2 as
below (I'm sure this can be written more cleanly, but you'll get the idea)
LpSolve2[a_,q_] := Block[{mat, newmat, newq=-q, len},
mat = Outer[Coefficient,
Flatten[p.a + Transpose[a] . p], Union@@p];
len = Sqrt[Length[mat]];
newmat = Table[Table[mat[[j]], {j, k*len+1, (k+1)*len}], {k, 0,
len-1}];
p/.Thread[(Union@@p) ->
LinearSolve[Flatten[Table[newmat[[j, k]], {j, 1, len}, {k, 1, j}],
1],
Flatten[Table[newq[[j, k]], {j, 1, Length[q]}, {k, 1, j}]]
]
] ]
then the precision of the result is much better. There is a corresponding
price paid in speed, though (for the bignum version).
I should mention that a hefty part of the timing of LpSolve2 (for these
examples) is in extracting the coeff matrix rather than in actually
running LinearSolve. Thus one might try to speed this up by finding it
without explicitly calculating p.a + Transpose[a] . p and then extracting
coeffs w.r.t. the vector p.
>
> Allan Hayes
> hay at haystack.demon.co.uk
>
Daniel Lichtblau, WRI