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
141
142
143
144
145
146
147
148
149
150
|
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2 Final//EN">
<HTML>
<HEAD>
<TITLE>PINV Moore-Penrose Pseudoinverse
</TITLE>
</HEAD>
<BODY>
<H2>PINV Moore-Penrose Pseudoinverse
</H2>
<P>
Section: <A HREF=sec_array.html> Array Generation and Manipulations </A>
<H3>Usage</H3>
Calculates the Moore-Penrose pseudoinverse of a matrix.
The general syntax for its use is
<PRE>
y = pinv(A,tol)
</PRE>
<P>
or for a default specification of the tolerance <code>tol</code>,
<PRE>
y = pinv(A)
</PRE>
<P>
For any <code>m x n</code> matrix <code>A</code>, the Moore-Penrose pseudoinverse
is the unique <code>n x m</code> matrix <code>B</code> that satisfies the following
four conditions
<UL>
<LI> <code>A B A = A</code>
</LI>
<LI> <code>B A B = B</code>
</LI>
<LI> <code>(A B)' = A B</code>
</LI>
<LI> <code>(B A)' = B A</code>
</LI>
</UL>
Also, it is true that <code>B y</code> is the minimum norm, least squares
solution to <code>A x = y</code>. The Moore-Penrose pseudoinverse is computed
from the singular value decomposition of <code>A</code>, with singular values
smaller than <code>tol</code> being treated as zeros. If <code>tol</code> is not specified
then it is chosen as
<PRE>
tol = max(size(A)) * norm(A) * teps(A).
</PRE>
<P>
<H3>Function Internals</H3>
The calculation of the MP pseudo-inverse is almost trivial once the
svd of the matrix is available. First, for a real, diagonal matrix
with positive entries, the pseudo-inverse is simply
<P>
<DIV ALIGN="CENTER">
<IMG SRC="pinv_eqn1.png">
</DIV>
<P>
One can quickly verify that this choice of matrix satisfies the
four properties of the pseudoinverse. Then, the pseudoinverse
of a general matrix <code>A = U S V'</code> is defined as
<P>
<DIV ALIGN="CENTER">
<IMG SRC="pinv_eqn2.png">
</DIV>
<P>
and again, using the facts that <code>U' U = I</code> and <code>V V' = I</code>, one
can quickly verify that this choice of pseudoinverse satisfies the
four defining properties of the MP pseudoinverse. Note that in
practice, the diagonal pseudoinverse <code>S^{+}</code> is computed with
a threshold (the <code>tol</code> argument to <code>pinv</code>) so that singular
values smaller than <code>tol</code> are treated like zeros.
<H3>Examples</H3>
Consider a simple <code>1 x 2</code> matrix example, and note the various
Moore-Penrose conditions:
<PRE>
--> A = float(rand(1,2))
A =
0.9526 0.4847
--> B = pinv(A)
B =
0.8338
0.4243
--> A*B*A
ans =
0.9526 0.4847
--> B*A*B
ans =
0.8338
0.4243
--> A*B
ans =
1.0000
--> B*A
ans =
0.7943 0.4042
0.4042 0.2057
</PRE>
<P>
To demonstrate that <code>pinv</code> returns the least squares solution,
consider the following very simple case
<PRE>
--> A = float([1;1;1;1])
A =
1
1
1
1
</PRE>
<P>
The least squares solution to <code>A x = b</code> is just <code>x = mean(b)</code>,
and computing the <code>pinv</code> of <code>A</code> demonstrates this
<PRE>
--> pinv(A)
ans =
0.2500 0.2500 0.2500 0.2500
</PRE>
<P>
Similarly, we can demonstrate the minimum norm solution with
the following simple case
<PRE>
--> A = float([1,1])
A =
1 1
</PRE>
<P>
The solutions of <code>A x = 5</code> are those <code>x_1</code> and <code>x_2</code> such that
<code>x_1 + x_2 = 5</code>. The norm of <code>x</code> is <code>x_1^ + x_2^2</code>, which is
<code>x_1^2 + (5-x_1)^2</code>, which is minimized for <code>x_1 = x_2 = 2.5</code>:
<PRE>
--> pinv(A) * 5.0
ans =
2.5000
2.5000
</PRE>
<P>
</BODY>
</HTML>
|
