-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathvandermonde.html
319 lines (286 loc) · 16.4 KB
/
vandermonde.html
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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
<!DOCTYPE html>
<html lang="en">
<head>
<!-- 2019-09-07 Sat 15:56 -->
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>Vandermonde matrices</title>
<meta name="generator" content="Org mode">
<meta name="author" content="George Kontsevich">
<meta name="description" content="Solving polynomials using Vandermonde matrices"
>
<link rel="stylesheet" type="text/css" href="../web/worg.css" />
<link rel="shortcut icon" href="../web/panda.svg" type="image/x-icon">
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
displayAlign: "center",
displayIndent: "0em",
"HTML-CSS": { scale: 100,
linebreaks: { automatic: "false" },
webFont: "TeX"
},
SVG: {scale: 100,
linebreaks: { automatic: "false" },
font: "TeX"},
NativeMML: {scale: 100},
TeX: { equationNumbers: {autoNumber: "AMS"},
MultLineWidth: "85%",
TagSide: "right",
TagIndent: ".8em"
}
});
</script>
<script type="text/javascript"
src="../MathJax/MathJax.js?config=TeX-AMS_CHTML"></script>
</head>
<body>
<div id="org-div-home-and-up">
<a accesskey="h" href="./index.html"> UP </a>
|
<a accesskey="H" href=".."> HOME </a>
</div><div id="content">
<h1 class="title">Vandermonde matrices</h1>
<div id="table-of-contents">
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#org9bdda0d">Solving polynomials</a></li>
<li><a href="#org54c0e08">Vandermonde matrices</a></li>
<li><a href="#org0d0c5da">Full rank case</a></li>
<li><a href="#org5d3367f">Least Squares</a></li>
</ul>
</div>
</div>
<div id="outline-container-org9bdda0d" class="outline-2">
<h2 id="org9bdda0d">Solving polynomials</h2>
<div class="outline-text-2" id="text-org9bdda0d">
<p>
<code>pages 185-186</code> introduces the Vandermonde matix which provides us with a way to solve for polynomials that fit a set of given points.
</p>
<p>
Polynomials are equations of the form:
</p>
\begin{equation}
y=a_{1}+a_{2}x+a_{3}x^{2}+a_{4}x^{3}+...
\end{equation}
<p>
It's a typicaly non-linear function of the form <b>y=f(x)</b> that takes inputs <code>x</code> and returns outputs <code>y</code>. Given a polynomial (ie. given its <b>a<sub>n</sub></b> factors) and given a value for <code>y</code> there can be multiple solutions for <code>x</code>. At first blush this doesn't seem related to the problems we've been tackling with matrices.
</p>
<p>
The Vandermonde matrix does a clever inversion of our typical problem. While we typically are looking for <code>x</code> given <code>y</code> here because we are ultimately interested in finding a polynomial <i>that fits</i> what we do instead is set up the problem to solve for the <code>a_{n}</code> polynomial factors. Instead of taking <b>y=f(x)</b> and solving for <b>x</b> we are taking a corresponding <b>y=f(a)</b> and solving for the polynomial factors <b>a<sub>n</sub></b>. We will see that the <b>x</b>'s are baked into the <b>f()</b>.
</p>
<p>
To see how this is built up we simply take the polynomial function and a series of points/measurements <code>[ (x_1,y_1) (x_2,y_2) (x_3,y_3) ... ]</code> and write out the equations.
</p>
\begin{equation}
y_1=a_{1}+a_{2}x_1+a_{3}x_{1}^{2}+a_{4}x_{1}^{3}+...\\
y_2=a_{1}+a_{2}x_2+a_{3}x_{2}^{2}+a_{4}x_{2}^{3}+...\\
y_3=a_{1}+a_{2}x_3+a_{3}x_{3}^{2}+a_{4}x_{3}^{3}+...\\
...
\end{equation}
<p>
No we we want to pull out the <b>a</b>'s to make this look like an <b>Ax=b</b> type of equation. While the <code>x</code>'s in the equations look nonlinear they're actually known values coordinates of our <code>x,y</code> points so they're fixed scalar values and they not something we are solving for. Writen out in matrix form we get
</p>
\begin{equation}
\begin{bmatrix}
1 & x_1 & x_{1}^2 & x_{1}^3 ..\\
1 & x_2 & x_{2}^2 & x_{2}^3 ..\\
1 & x_3 & x_{3}^2 & x_{3}^3 ..\\
...\\
\end{bmatrix}
\begin{bmatrix}
a_1\\
a_2\\
a_3\\
a_4\\
...\\
\end{bmatrix}
=
\begin{bmatrix}
y_1\\
y_2\\
y_3\\
...\\
\end{bmatrix}
\end{equation}
<p>
Again, all the values of the x-matrix are known. (the x-matrix is the actual <i>Vandermonde matrix</i> - <b>V</b>)
</p>
</div>
</div>
<div id="outline-container-org54c0e08" class="outline-2">
<h2 id="org54c0e08">Vandermonde matrices</h2>
<div class="outline-text-2" id="text-org54c0e08">
<div class="org-src-container">
<pre class="src src-clojure">(<span style="color: #0000FF;">ns</span> <span style="color: #6434A3;">morelinear.vandermonde</span>
(<span style="color: #D0372D;">:require</span> [morelinear.leastsquares <span style="color: #D0372D;">:as</span> leastsquares<span style="color: #999999;">]</span>
[clojure.core.matrix <span style="color: #D0372D;">:refer</span> <span style="color: #D0372D;">:all</span><span style="color: #999999;">]</span>
[clojure.core.matrix.linear <span style="color: #D0372D;">:as</span> linear<span style="color: #999999;">]))</span>
(set-current-implementation <span style="color: #D0372D;">:vectorz</span>)
</pre>
</div>
<p>
So given a series of <code>x</code>'s we can build up this matrix directly.
</p>
<div class="org-src-container">
<pre class="src src-clojure">(<span style="color: #0000FF;">defn</span> <span style="color: #006699;">polynomial-matrix</span>
<span style="color: #036A07;">"Take a vector of input X's and build a vandermonde matrix.</span>
<span style="color: #036A07;"> The DEGREE shouldn't be larger than the number of values given.</span>
<span style="color: #036A07;"> If non provided it will be chosen to be the maximal value</span>
<span style="color: #036A07;"> and your resulting matrix will be square"</span>
([polynomial-input-values
degree<span style="color: #999999;">]</span>
(mapv (<span style="color: #0000FF;">fn</span> [x<span style="color: #999999;">]</span>
(mapv (<span style="color: #0000FF;">fn</span> [power<span style="color: #999999;">]</span>
(pow x power<span style="color: #999999;">))</span>
(range degree<span style="color: #999999;">)))</span>
polynomial-input-values<span style="color: #999999;">))</span>
<span style="color: #8D8D84;">;; </span><span style="color: #8D8D84; font-style: italic;">Default square matrix case</span>
([polynomial-input-values<span style="color: #999999;">]</span>
(polynomial-matrix polynomial-input-values
(ecount polynomial-input-values<span style="color: #999999;">))))</span>
</pre>
</div>
<p>
The number of polynomials terms (ie. the <b>degree</b> <i>n</i> in the last x<sup>n</sup>) is variable.
</p>
</div>
</div>
<div id="outline-container-org0d0c5da" class="outline-2">
<h2 id="org0d0c5da">Full rank case</h2>
<div class="outline-text-2" id="text-org0d0c5da">
<p>
If we have as many polynomial terms as we have points then we are left with a square x-matrix which we can solve using Gaussian elimination (excluding singular/degenerate cases). In other words we can use the <b>LU decomposition</b> to solve for all the polynomial factors <b>a<sub>n</sub></b> in this equivalent system <b>Va=y</b> - where the <b>V</b> is the matrix of exponents of <code>x</code>.
</p>
<p>
Once we solve for the <b>a<sub>n</sub></b>'s we can then stick them back into non-linear polynomial equation we started looking at in the beginning and get a solution.
</p>
\begin{equation}
y=a_{1}+a_{2}x+a_{3}x^{2}+a_{4}x^{3}+...
\end{equation}
<p>
It will not only hold true for all our <code>(x,y)</code> points but also for all other values of <code>x</code> we want to test - so we will have in the end fit a polynomial curve through all our points.
</p>
</div>
</div>
<div id="outline-container-org5d3367f" class="outline-2">
<h2 id="org5d3367f">Least Squares</h2>
<div class="outline-text-2" id="text-org5d3367f">
<p>
However with our least squares method we now have an alternate solution. We can choose a <b>degree</b> that is <i>lower</i> than the number of points/values we have and then solve for approximate solutions. This is often the preferrable solution. In general when we fit a polynomial we have a decent idea of the degree of the function that drives whatever we are measuring. We also know that the input values are not perfect, so if we overfit the points with a high-degree polynomial the results will also be very noisy. So typically we preselect the degree of the polynomial and then try to give the function as many data points as we can to get the best approximation of the underlying function as we can.
</p>
<blockquote>
<p>
<b>Note</b>: That this doesn't always hold. If you have some non-linear function that you are approximating with a Taylor series, then cut off the number of polynomials serves no purpose
</p>
</blockquote>
<p>
As before the Least Squares solutions come in a few flavors and it also works on the square case (even when the matrix is singular) - so we can skip writing a function for the square case and just have one polynomial solver for all square and skinny matrices
</p>
<div class="org-src-container">
<pre class="src src-clojure">(<span style="color: #0000FF;">defn</span> <span style="color: #006699;">split-into-x-y</span>
[points<span style="color: #999999;">]</span>
(<span style="color: #0000FF;">let</span> [vector-of-xs-and-ys (apply mapv vector points<span style="color: #999999;">)]</span>
{<span style="color: #D0372D;">:x</span> (mutable (first vector-of-xs-and-ys<span style="color: #999999;">))</span>
<span style="color: #D0372D;">:y</span> (mutable (second vector-of-xs-and-ys<span style="color: #999999;">))}))</span>
(<span style="color: #0000FF;">defmulti</span> <span style="color: #006699;">polynomial-factors</span>
<span style="color: #036A07;">"Solve for a polynomial equation of DEGREE that fits the given POINTS</span>
<span style="color: #036A07;"> and return the list of polynomial factors.</span>
<span style="color: #036A07;"> the underlying least squares METHOD used can be:</span>
<span style="color: #036A07;"> :lu</span>
<span style="color: #036A07;"> :lu-jumbo</span>
<span style="color: #036A07;"> :householder"</span>
(<span style="color: #0000FF;">fn</span> [degree
points
method<span style="color: #999999;">]</span>
method<span style="color: #999999;">))</span>
(<span style="color: #0000FF;">defmethod</span> <span style="color: #006699;">polynomial-factors</span> <span style="color: #D0372D;">:lu</span>
[degree
points
method<span style="color: #999999;">]</span>
(<span style="color: #0000FF;">let</span> [xy (split-into-x-y points<span style="color: #999999;">)]</span>
(<span style="color: #6434A3;">leastsquares</span>/lu-direct (polynomial-matrix (<span style="color: #D0372D;">:x</span> xy<span style="color: #999999;">)</span>
degree<span style="color: #999999;">)</span>
(<span style="color: #D0372D;">:y</span> xy<span style="color: #999999;">))))</span>
(<span style="color: #0000FF;">defmethod</span> <span style="color: #006699;">polynomial-factors</span> <span style="color: #D0372D;">:lu-jumbo</span>
[degree
points
method<span style="color: #999999;">]</span>
(<span style="color: #0000FF;">let</span> [xy (split-into-x-y points<span style="color: #999999;">)]</span>
(<span style="color: #6434A3;">leastsquares</span>/lu-jumbo (polynomial-matrix (<span style="color: #D0372D;">:x</span> xy<span style="color: #999999;">)</span>
degree<span style="color: #999999;">)</span>
(<span style="color: #D0372D;">:y</span> xy<span style="color: #999999;">))))</span>
(<span style="color: #0000FF;">defmethod</span> <span style="color: #006699;">polynomial-factors</span> <span style="color: #D0372D;">:householder</span>
[degree
points
method<span style="color: #999999;">]</span>
(<span style="color: #0000FF;">let</span> [xy (split-into-x-y points<span style="color: #999999;">)</span>
vandermonde-matrix (mutable (polynomial-matrix (<span style="color: #D0372D;">:x</span> xy<span style="color: #999999;">)</span>
degree<span style="color: #999999;">))]</span>
(<span style="color: #6434A3;">leastsquares</span>/householder-qr vandermonde-matrix
(<span style="color: #D0372D;">:y</span> xy<span style="color: #999999;">))))</span>
</pre>
</div>
<div class="org-src-container">
<pre class="src src-clojure">(polynomial-factors 3
[[138.0 165.0<span style="color: #999999;">]</span>
[275.0 197.0<span style="color: #999999;">]</span>
[277.0 262.0<span style="color: #999999;">]</span>
[186.0 277.0<span style="color: #999999;">]</span>
[154.0 240.0<span style="color: #999999;">]</span>
[236.0 141.0<span style="color: #999999;">]</span>
[322.0 123.0<span style="color: #999999;">]</span>
[358.0 122.0<span style="color: #999999;">]</span>
[182.0 289.0<span style="color: #999999;">]</span>
[118.0 294.0<span style="color: #999999;">]</span>
[108.0 292.0<span style="color: #999999;">]</span>
[54.0 259.0<span style="color: #999999;">]]</span>
<span style="color: #D0372D;">:householder</span><span style="color: #999999;">)</span>
<span style="color: #8D8D84;">;;</span><span style="color: #8D8D84; font-style: italic;">(251.1770462843548 0.3111353706549809 -0.001904254283169026)</span>
</pre>
</div>
<p>
Now that we've got the polynomial factors we want to just wrap this up and be able to get the polynomial function itself
</p>
<div class="org-src-container">
<pre class="src src-clojure">(<span style="color: #0000FF;">defn</span> <span style="color: #006699;">polynomial-function</span>
[degree points method<span style="color: #999999;">]</span>
(<span style="color: #0000FF;">let</span> [factors (polynomial-factors degree
points
method<span style="color: #999999;">)]</span>
(<span style="color: #0000FF;">fn</span> [x] (reduce + (map-indexed (<span style="color: #0000FF;">fn</span> [index
factor<span style="color: #999999;">]</span>
(* factor
(pow x index<span style="color: #999999;">)))</span>
factors<span style="color: #999999;">)))))</span>
</pre>
</div>
<div class="org-src-container">
<pre class="src src-clojure">(<span style="color: #0000FF;">let</span> [our-function (polynomial-function 3
[[138.0 165.0<span style="color: #999999;">]</span>
[275.0 197.0<span style="color: #999999;">]</span>
[277.0 262.0<span style="color: #999999;">]</span>
[186.0 277.0<span style="color: #999999;">]</span>
[154.0 240.0<span style="color: #999999;">]</span>
[236.0 141.0<span style="color: #999999;">]</span>
[322.0 123.0<span style="color: #999999;">]</span>
[358.0 122.0<span style="color: #999999;">]</span>
[182.0 289.0<span style="color: #999999;">]</span>
[118.0 294.0<span style="color: #999999;">]</span>
[108.0 292.0<span style="color: #999999;">]</span>
[54.0 259.0<span style="color: #999999;">]]</span>
<span style="color: #D0372D;">:householder</span><span style="color: #999999;">)]</span>
(println (our-function 11<span style="color: #999999;">))</span>
(println (our-function 50<span style="color: #999999;">))</span>
(println (our-function 70<span style="color: #999999;">))</span>
(println (our-function 150<span style="color: #999999;">))</span>
(println (our-function 250<span style="color: #999999;">))</span>
(println (our-function 350<span style="color: #999999;">)))</span>
<span style="color: #8D8D84;">;;</span><span style="color: #8D8D84; font-style: italic;">(251.1770462843548 0.3111353706549809 -0.001904254283169026)</span>
</pre>
</div>
</div>
</div>
</div>
</body>
</html>