[Date Index]
[Thread Index]
[Author Index]
Re: RowReduce of badly conditioned matrix
*To*: mathgroup at smc.vnet.net
*Subject*: [mg25248] Re: [mg25227] RowReduce of badly conditioned matrix
*From*: Daniel Lichtblau <danl at wolfram.com>
*Date*: Sun, 17 Sep 2000 04:47:27 -0400 (EDT)
*References*: <200009150621.CAA10220@smc.vnet.net>
*Sender*: owner-wri-mathgroup at wolfram.com
Camillo Ressl wrote:
>
> Hallo!
>
> I have the following problem: I have a matrix a with 24 rows and 42
> colums.Its elements are given in double precision. For this matrix the
> following holds:
>
> Rank[a] == Rank[a.Transpose[a]] = 21
>
> In[298]:= {u,md,v}= SingularValues[a];
>
> In[299]:= md
>
> Out[299]= {1.41421, 1.41421, 1.41421, 1.13491, 1.13488, 1.00081,
> 1.00002,
>
> > 1.00001, 1.00001, 1.00001, 1., 1., 1., 0.999999, 0.999987,
> 0.999987,
>
> > 0.999986, 0.999978, 0.998748, 0.844348, 0.843791}
>
> And now I want to compute the rowreduced matrix of a and I get the
> following error message.
>
> In[300]:= RowReduce[a]
>
> RowReduce::luc:
> Result for RowReduce of badly conditioned matrix {<<24>>}
> may contain significant numerical errors.
>
> 7
> Out[300]= {{1, 0, 0, 0, 0, 0, 0, 0, 4.78983 10 , 0, 0, 0, 0, 0, 0, 0,
> ...etc.
>
> I have no idea what the reason for this message could be since the
> Singular Values seem o.k.
> And the rank-defect of 3 should not be the reason since I expect 3
> zero-rows for the rowreduced matrix.
>
> Is there someone who can help me or give me at least a few hints?
>
> Many thanks for any help.
>
> Camillo
The problem is that RowReduce used a "small" but nonzero pivot. The
warning is quite useful, I think, since it alerts you to the fact that
something is possibly amiss.
The reason for this behavior is simple. RowReduce uses Gaussian
elimination with partial pivoting, and this has essentially no ability
to cope with numerical problems associated with rank.
One possible workaround is to provide a nondefault ZeroTest to
RowReduce. This is problematic because it is not impervious to scaling
of the input. That said, here is a quick demo along the lines of the
example you outline.
We generate a 21x42 matrix, then generate three random linear
combinations of its rows and join those to the matrix to form a 24x42
matrix of rank 21.
In[44]:= mat = Table[Random[], {21}, {42}];
In[45]:= extras = Table[Apply[Plus,Table[Random[],{21}]*mat], {3}];
In[46]:= newmat = Join[mat, extras];
We use SingularValues to check that rank.
In[47]:= {uu,md,vv}= SingularValues[newmat];
In[48]:= Length[md]
Out[48]= 21
Against our better judgement, we'll try RowReduce.
In[52]:= dontdothis = RowReduce[newmat];
RowReduce::luc:
Result for RowReduce of badly conditioned matrix {<<24>>}
may contain significant numerical errors.
In[53]:= Max[Abs[Take[dontdothis,-3]]]
Out[53]= 12.7634
Here is how one might row reduce using the ZeroTest. We check that the
last three rows are all zeroes.
In[56]:= rred = RowReduce[newmat, ZeroTest->(Chop[#]==0&)];
In[57]:= Max[Abs[Take[rred,-3]]]
Out[57]= 0
As I mentioned, this is still not a really good way to proceed. A more
numerically sound approach is to preprocess with a QR decomposition, as
this tends to discover deficiencies in rank.
In[59]:= {qq,rr} = QRDecomposition[newmat];
Now row reduce the upper triangular part (the "R" is "QR").
In[60]:= rred2 = RowReduce[rr];
We'll compare to the first 21 rows of rred. We get results that agree to
within 2-3 orders of magnitude of machine epsilon.
In[61]:= Max[Abs[Drop[rred,-3] - rred2]]
-14
Out[61]= 8.26006 10
This QR method is useful if you do not want to go to the trouble of a
singular values decomposition. In the case where you have in fact done
that (as we did above), you can work with the svd instead. In essence we
just discard the premultiplier matrix uu.
In[68]:= rred3 = RowReduce[md*vv];
In[69]:= Max[Abs[rred3 - rred2]]
-14
Out[69]= 1.59872 10
Daniel Lichtblau
Wolfram Research
Prev by Date:
**Re: Displaying Mixed Numbers**
Next by Date:
**point inside torus**
Previous by thread:
**RowReduce of badly conditioned matrix**
Next by thread:
**ParametricPlot3D is buggy**
| |