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
|
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2 Final//EN">
<HTML>
<HEAD>
<TITLE>LEFTDIVIDE Matrix Equation Solver/Divide Operator
</TITLE>
</HEAD>
<BODY>
<H2>LEFTDIVIDE Matrix Equation Solver/Divide Operator
</H2>
<P>
Section: <A HREF=sec_operators.html> Mathematical Operators </A>
<H3>Usage</H3>
The divide operator <code>\</code> is really a combination of three
operators, all of which have the same general syntax:
<PRE>
Y = A \ B
</PRE>
<P>
where <code>A</code> and <code>B</code> are arrays of numerical type. The result <code>Y</code> depends
on which of the following three situations applies to the arguments
<code>A</code> and <code>B</code>:
<OL>
<LI> <code>A</code> is a scalar, <code>B</code> is an arbitrary <code>n</code>-dimensional numerical array, in which case the output is each element of <code>B</code> divided by the scalar <code>A</code>.
</LI>
<LI> <code>A,B</code> are matrices with the same number of rows, i.e., <code>A</code> is of size <code>M x K</code>, and <code>B</code> is of size <code>M x L</code>, in which case the output is of size <code>K x L</code>.
</LI>
</OL>
The output follows the standard type promotion rules, although in the first two cases, if <code>A</code> and <code>B</code> are integers, the output is an integer also, while in the third case if <code>A</code> and <code>B</code> are integers, the output is of type <code>double</code>.
A few additional words about the third version, in which <code>A</code> and <code>B</code> are matrices. Very loosely speaking, <code>Y</code> is the matrix that satisfies <code>A * Y = B</code>. In cases where such a matrix exists. If such a matrix does not exist, then a matrix <code>Y</code> is returned that approximates <code>A * Y \approx B</code>.
<H3>Function Internals</H3>
There are three formulae for the times operator. For the first form
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn1.png">
</DIV>
<P>
In the second form, the calculation of the output depends on the size of <code>A</code>. Because each column of <code>B</code> is treated independantly, we can rewrite the equation <code>A Y = B</code> as
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn2.png">
</DIV>
<P>
where <code>y_i</code> are the columns of <code>Y</code>, and <code>b_i</code> are the columns of the matrix <code>B</code>. If <code>A</code> is a square matrix, then the LAPACK routine <code>*gesvx</code> (where the <code>*</code> is replaced with <code>sdcz</code> depending on the type of the arguments) is used, which uses an LU decomposition of <code>A</code> to solve the sequence of equations sequentially. If <code>A</code> is singular, then a warning is emitted.
On the other hand, if <code>A</code> is rectangular, then the LAPACK routine <code>*gelsy</code> is used. Note that these routines are designed to work with matrices <code>A</code> that are full rank - either full column rank or full row rank. If <code>A</code> fails to satisfy this assumption, a warning is emitted. If <code>A</code> has full column rank (and thus necessarily has more rows than columns), then theoretically, this operator finds the columns <code>y_i</code> that satisfy:
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn3.png">
</DIV>
<P>
and each column is thus the Least Squares solution of <code>A y = b_i</code>. On the other hand, if <code>A</code> has full row rank (and thus necessarily has more columns than rows), then theoretically, this operator finds the columns <code>y_i</code> that satisfy
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn4.png">
</DIV>
<P>
and each column is thus the Minimum Norm vector <code>y_i</code> that satisfies <code>A y_i = b_i</code>.
In the event that the matrix <code>A</code> is neither full row rank nor full column rank, a solution is returned, that is the minimum norm least squares solution. The solution is computed using an orthogonal factorization technique that is documented in the LAPACK User's Guide (see the References section for details).
<H3>Examples</H3>
Here are some simple examples of the divide operator. We start with a simple example of a full rank, square matrix:
<PRE>
--> A = [1,1;0,1]
A =
1 1
0 1
</PRE>
<P>
Suppose we wish to solve
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn5.png">
</DIV>
<P>
(which by inspection has the solution <code>y_1 = 1</code>, <code>y_2 = 2</code>). Thus we compute:
<PRE>
--> B = [3;2]
B =
3
2
--> Y = A\B
Y =
1
2
</PRE>
<P>
Suppose we wish to solve a trivial Least Squares (LS) problem. We want to find a simple scaling of the vector <code>[1;1]</code> that is closest to the point <code>[2,1]</code>. This is equivalent to solving
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn6.png">
</DIV>
<P>
in a least squares sense. For fun, we can calculate the solution using calculus by hand. The error we wish to minimize is
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn7.png">
</DIV>
<P>
Taking a derivative with respect to <code>y</code>, and setting to zero (which we must have for an extrema when <code>y</code> is unconstrained)
<P>
<DIV ALIGN="CENTER">
<IMG SRC="leftdivide_eqn8.png">
</DIV>
<P>
which we can simplify to <code>4y = 6</code> or <code>y = 3/2</code> (we must, technically, check to make sure this is a minimum, and not a maximum or an inflection point). Here is the same calculation performed using FreeMat:
<PRE>
--> A = [1;1]
A =
1
1
--> B = [2;1]
B =
2
1
--> A\B
ans =
1.5000
</PRE>
<P>
which is the same solution.
</BODY>
</HTML>
|
