1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
<!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN">
<!--Converted with LaTeX2HTML 96.1-h (September 30, 1996) by Nikos Drakos (nikos@cbl.leeds.ac.uk), CBLU, University of Leeds -->
<HTML>
<HEAD>
<TITLE>Error Bounds for Linear Least Squares Problems</TITLE>
<META NAME="description" CONTENT="Error Bounds for Linear Least Squares Problems">
<META NAME="keywords" CONTENT="slug">
<META NAME="resource-type" CONTENT="document">
<META NAME="distribution" CONTENT="global">
<LINK REL=STYLESHEET HREF="slug.css">
</HEAD>
<BODY LANG="EN" >
<A NAME="tex2html3962" HREF="node141.html"><IMG WIDTH=37 HEIGHT=24 ALIGN=BOTTOM ALT="next" SRC="http://www.netlib.org/utk/icons/next_motif.gif"></A> <A NAME="tex2html3960" HREF="node132.html"><IMG WIDTH=26 HEIGHT=24 ALIGN=BOTTOM ALT="up" SRC="http://www.netlib.org/utk/icons/up_motif.gif"></A> <A NAME="tex2html3954" HREF="node139.html"><IMG WIDTH=63 HEIGHT=24 ALIGN=BOTTOM ALT="previous" SRC="http://www.netlib.org/utk/icons/previous_motif.gif"></A> <A NAME="tex2html3964" HREF="node1.html"><IMG WIDTH=65 HEIGHT=24 ALIGN=BOTTOM ALT="contents" SRC="http://www.netlib.org/utk/icons/contents_motif.gif"></A> <A NAME="tex2html3965" HREF="node190.html"><IMG WIDTH=43 HEIGHT=24 ALIGN=BOTTOM ALT="index" SRC="http://www.netlib.org/utk/icons/index_motif.gif"></A> <BR>
<B> Next:</B> <A NAME="tex2html3963" HREF="node141.html">Error Bounds for the </A>
<B>Up:</B> <A NAME="tex2html3961" HREF="node132.html">Accuracy and Stability</A>
<B> Previous:</B> <A NAME="tex2html3955" HREF="node139.html">Error Bounds for Linear </A>
<BR> <P>
<H1><A NAME="SECTION04660000000000000000">Error Bounds for Linear Least Squares Problems</A></H1>
<A NAME="seclsq"> </A>
<P>
The linear least squares problem is to find <I>x</I> that minimizes
<IMG WIDTH=73 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18475" SRC="img610.gif">.
We discuss error bounds for the most common case where <I>A</I> is <I>m</I>-by-<I>n</I>
with <I>m</I> > <I>n</I>, and <I>A</I> has full rank<A NAME="5322"> </A>;
this is called an <EM>overdetermined least squares problem</EM>
<A NAME="5324"> </A>
(the following code fragments deal with <I>m</I>=<I>n</I> as well).
<P>
Let <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> be the solution computed by the driver routine
PxGELS (see section <A HREF="node45.html#subsecdrivellsq">3.2.2</A>).
An approximate error
bound
<A NAME="5328"> </A><A NAME="5329"> </A><A NAME="5330"> </A><A NAME="5331"> </A>
<BR><IMG WIDTH=318 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18463" SRC="img611.gif"><BR>
may be computed in the following way:
<P>
<PRE> EPSMCH = PSLAMCH( ICTXT, 'E' )
* Get the 2-norm of the right hand side B
BNORM = PSLANGE( 'F', M, 1, B, IB, JB, DESCB, WORK )
* Solve the least squares problem; the solution X overwrites B
CALL PSGELS( 'N', M, N, 1, A, IA, JA, DESCA, B, IB, JB, DESCB,
$ WORK, LWORK, INFO )
IF( MIN( M, N ).GT.0 ) THEN
* Get the 2-norm of the residual A*X-B
RNORM = PSLANGE( 'F', M-N, 1, B, IB+N, JB, DESCB, WORK )
* Get the reciprocal condition number RCOND of A
CALL PSTRCON( 'I', 'U', 'N', N, A, IA, JA, DESCA, RCOND, WORK,
$ LWORK, IWORK, LIWORK, INFO )
RCOND = MAX( RCOND, EPSMCH )
IF( BNORM.GT.0.0 ) THEN
SINT = RNORM / BNORM
ELSE
SINT = 0.0
END IF
COST = MAX( SQRT( ( 1.0E0 - SINT )*( 1.0E0 + SINT ) ), EPSMCH )
TANT = SINT / COST
ERRBD = EPSMCH*( 2.0E0 / ( RCOND*COST ) + TANT / RCOND**2 )
END IF</PRE>
<A NAME="5335"> </A>
<P>
For example, if <IMG WIDTH=315 HEIGHT=30 ALIGN=MIDDLE ALT="tex2html_wrap_inline18242" SRC="img571.gif">,
<BR><IMG WIDTH=406 HEIGHT=87 ALIGN=BOTTOM ALT="displaymath18464" SRC="img612.gif"><BR>
then, to four decimal places,
<BR><IMG WIDTH=340 HEIGHT=67 ALIGN=BOTTOM ALT="displaymath18465" SRC="img613.gif"><BR>
<IMG WIDTH=106 HEIGHT=13 ALIGN=BOTTOM ALT="tex2html_wrap_inline18495" SRC="img614.gif">,
<IMG WIDTH=107 HEIGHT=13 ALIGN=BOTTOM ALT="tex2html_wrap_inline18497" SRC="img615.gif">,
<IMG WIDTH=155 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18499" SRC="img616.gif">,
<IMG WIDTH=137 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18501" SRC="img617.gif">, and the true error
is <IMG WIDTH=68 HEIGHT=16 ALIGN=BOTTOM ALT="tex2html_wrap_inline18503" SRC="img618.gif">.
<P>
Note that in the preceding code fragment, the routine PSLANGE was
used to compute the two-norm of the right hand side matrix <I>B</I> and
the residual <IMG WIDTH=79 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline18507" SRC="img619.gif">. This routine was chosen because the result
of the computation (BNORM or RNORM, respectively) is automatically
known on all process columns within the process grid. The routine
PSNRM2 could have also been
used to perform this calculation; however, the use of PSNRM2 in this
example would have required an additional communication broadcast,
because the resulting value of BNORM or RNORM, respectively, is known only
within the process column owning B(:,JB).
<P>
<B>Further Details: Error Bounds for Linear Least Squares
Problems</B><A NAME="secbackgroundlsq"> </A>
<P>
The conventional error analysis of linear least squares problems goes
as follows<A NAME="5353"> </A>.
As above, let <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> be the solution to minimizing
<IMG WIDTH=73 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18475" SRC="img610.gif"> computed by
ScaLAPACK using the least squares driver PxGELS
(see section <A HREF="node45.html#subsecdrivellsq">3.2.2</A>).
We discuss the most common case, where <I>A</I> is
overdetermined<A NAME="5355"> </A>
(i.e., has more rows than columns) and has full rank [<A HREF="node189.html#GVL2">71</A>]:
<A NAME="5357"> </A><A NAME="5358"> </A><A NAME="5359"> </A><A NAME="5360"> </A>
<P>
<BLOCKQUOTE> The computed solution <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> has a small normwise backward error.
In other words, <IMG WIDTH=9 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17482" SRC="img449.gif"> minimizes <IMG WIDTH=167 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18519" SRC="img620.gif">, where
<I>E</I> and <I>f</I> satisfy
<A NAME="5362"> </A>
<A NAME="5363"> </A>
<BR><IMG WIDTH=355 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18466" SRC="img621.gif"><BR>
and <I>p</I>(<I>n</I>) is a modestly growing function of <I>n</I>. We take <I>p</I>(<I>n</I>)=1 in
the code fragments above.
Let <IMG WIDTH=195 HEIGHT=27 ALIGN=MIDDLE ALT="tex2html_wrap_inline18531" SRC="img622.gif"> (approximated by
1/<TT>RCOND</TT> in the above code fragments),
<IMG WIDTH=107 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18533" SRC="img623.gif"> (= <TT>RNORM</TT> above), and <IMG WIDTH=115 HEIGHT=27 ALIGN=MIDDLE ALT="tex2html_wrap_inline18535" SRC="img624.gif">
(<TT>SINT = RNORM / BNORM</TT> above). Here, <IMG WIDTH=7 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline17638" SRC="img485.gif"> is the acute angle between
the vectors <IMG WIDTH=22 HEIGHT=12 ALIGN=BOTTOM ALT="tex2html_wrap_inline18539" SRC="img625.gif"> and <I>b</I>.
<A NAME="5375"> </A>
<A NAME="5376"> </A>
Then when <IMG WIDTH=41 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18543" SRC="img626.gif"> is small, the error <IMG WIDTH=40 HEIGHT=22 ALIGN=MIDDLE ALT="tex2html_wrap_inline18545" SRC="img627.gif"> is bounded by
<BR><IMG WIDTH=417 HEIGHT=41 ALIGN=BOTTOM ALT="displaymath18467" SRC="img628.gif"><BR>
where <IMG WIDTH=44 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18549" SRC="img629.gif"> = <TT>COST</TT> and <IMG WIDTH=45 HEIGHT=26 ALIGN=MIDDLE ALT="tex2html_wrap_inline18551" SRC="img630.gif"> = <TT>TANT</TT> in the code fragments
above.
</BLOCKQUOTE>
<P>
We avoid overflow by making sure <TT>RCOND</TT> and <TT>COST</TT> are both at least
<IMG WIDTH=24 HEIGHT=8 ALIGN=BOTTOM ALT="tex2html_wrap_inline18308" SRC="img594.gif"> <TT>EPSMCH</TT>, and by handling the case of a zero <TT>B</TT> matrix
separately (<TT>BNORM = 0</TT>).
<A NAME="5389"> </A>
<A NAME="5390"> </A>
<P>
<IMG WIDTH=195 HEIGHT=27 ALIGN=MIDDLE ALT="tex2html_wrap_inline18531" SRC="img622.gif"> may be computed directly
from the singular values of <I>A</I> returned by PxGESVD.
It may also be approximated by using PxTRCON following calls to
PxGELS. PxTRCON estimates <IMG WIDTH=22 HEIGHT=17 ALIGN=MIDDLE ALT="tex2html_wrap_inline18559" SRC="img631.gif"> or <IMG WIDTH=15 HEIGHT=17 ALIGN=MIDDLE ALT="tex2html_wrap_inline18561" SRC="img632.gif"> instead
of <IMG WIDTH=15 HEIGHT=17 ALIGN=MIDDLE ALT="tex2html_wrap_inline18563" SRC="img633.gif">, but these can differ from <IMG WIDTH=15 HEIGHT=17 ALIGN=MIDDLE ALT="tex2html_wrap_inline18563" SRC="img633.gif"> by at most a factor of <I>n</I>.
<A NAME="5394"> </A><A NAME="5395"> </A><A NAME="5396"> </A><A NAME="5397"> </A>
<A NAME="5398"> </A><A NAME="5399"> </A><A NAME="5400"> </A><A NAME="5401"> </A>
<P>
<HR><A NAME="tex2html3962" HREF="node141.html"><IMG WIDTH=37 HEIGHT=24 ALIGN=BOTTOM ALT="next" SRC="http://www.netlib.org/utk/icons/next_motif.gif"></A> <A NAME="tex2html3960" HREF="node132.html"><IMG WIDTH=26 HEIGHT=24 ALIGN=BOTTOM ALT="up" SRC="http://www.netlib.org/utk/icons/up_motif.gif"></A> <A NAME="tex2html3954" HREF="node139.html"><IMG WIDTH=63 HEIGHT=24 ALIGN=BOTTOM ALT="previous" SRC="http://www.netlib.org/utk/icons/previous_motif.gif"></A> <A NAME="tex2html3964" HREF="node1.html"><IMG WIDTH=65 HEIGHT=24 ALIGN=BOTTOM ALT="contents" SRC="http://www.netlib.org/utk/icons/contents_motif.gif"></A> <A NAME="tex2html3965" HREF="node190.html"><IMG WIDTH=43 HEIGHT=24 ALIGN=BOTTOM ALT="index" SRC="http://www.netlib.org/utk/icons/index_motif.gif"></A> <BR>
<B> Next:</B> <A NAME="tex2html3963" HREF="node141.html">Error Bounds for the </A>
<B>Up:</B> <A NAME="tex2html3961" HREF="node132.html">Accuracy and Stability</A>
<B> Previous:</B> <A NAME="tex2html3955" HREF="node139.html">Error Bounds for Linear </A>
<P><ADDRESS>
<I>Susan Blackford <BR>
Tue May 13 09:21:01 EDT 1997</I>
</ADDRESS>
</BODY>
</HTML>
|
