Skip to content

Commit 46d7cca

Browse files
fs446fs446
committed
line fit example of updated ex05
1 parent b648e9b commit 46d7cca

File tree

2 files changed

+341
-0
lines changed

2 files changed

+341
-0
lines changed

‎index.ipynb

+1
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@
4949
"## Exercise: Left Inverse / Linear Regression with Ordinary Least Squares (OLS)\n",
5050
"- [SVD and Left Inverse](exercise04_leftinv.ipynb)\n",
5151
"- [SVD and Right Inverse](exercise04_rightinv.ipynb)\n",
52+
"- [Line Fit with Linear Regression](line_fit_linear_regression.ipynb)\n",
5253
"- [Linear Regression with OLS](ols.ipynb)\n",
5354
"\n",
5455
"\n",

‎line_fit_linear_regression.ipynb

+340
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,340 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"Sascha Spors,\n",
8+
"Professorship Signal Theory and Digital Signal Processing,\n",
9+
"Institute of Communications Engineering (INT),\n",
10+
"Faculty of Computer Science and Electrical Engineering (IEF),\n",
11+
"University of Rostock,\n",
12+
"Germany\n",
13+
"\n",
14+
"# Data Driven Audio Signal Processing - A Tutorial with Computational Examples\n",
15+
"\n",
16+
"Winter Semester 2024/25 (Master Course #24512)\n",
17+
"\n",
18+
"- lecture: https://github.com/spatialaudio/data-driven-audio-signal-processing-lecture\n",
19+
"- tutorial: https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise\n",
20+
"\n",
21+
"Feel free to contact lecturer [email protected]"
22+
]
23+
},
24+
{
25+
"cell_type": "markdown",
26+
"metadata": {},
27+
"source": [
28+
"# Exercise 5: Linear Regression Toy Example\n",
29+
"\n",
30+
"## Objectives\n",
31+
"\n",
32+
"When no assumption on an underlying data generation process is being made, pure linear algebra is used to solve for model parameters. Hence, we should link\n",
33+
"- linear regression model (simple line fit)\n",
34+
"- left inverse of a tall / thin, full column (feature) matrix\n",
35+
"- (residual) least squares\n",
36+
"- projection matrices to the 4 subspaces\n",
37+
"\n",
38+
"to the very same playground using the following simple toy example."
39+
]
40+
},
41+
{
42+
"cell_type": "code",
43+
"execution_count": null,
44+
"metadata": {},
45+
"outputs": [],
46+
"source": [
47+
"import matplotlib.pyplot as plt\n",
48+
"import numpy as np\n",
49+
"from scipy.linalg import svd, diagsvd, inv, pinv, norm\n",
50+
"from numpy.linalg import matrix_rank"
51+
]
52+
},
53+
{
54+
"cell_type": "code",
55+
"execution_count": null,
56+
"metadata": {},
57+
"outputs": [],
58+
"source": [
59+
"X = np.array([[1, 1],\n",
60+
" [1, 2],\n",
61+
" [1, 3],\n",
62+
" [1, 4]])\n",
63+
"print(X, X.shape, matrix_rank(X))\n",
64+
"y_col = np.array([[1],\n",
65+
" [3],\n",
66+
" [5],\n",
67+
" [7]])\n",
68+
"print(y_col, y_col.shape)\n",
69+
"[U, s, Vh] = svd(X)\n",
70+
"V = Vh.T\n",
71+
"y_left_null = (-U[:,2]+U[:,3])[:, None] # [:, None] makes it a (4,1) array\n",
72+
"print(y_left_null, y_left_null.shape)\n",
73+
"y = y_col + y_left_null\n",
74+
"print(y, y.shape)\n",
75+
"M, N = X.shape\n",
76+
"print(M, N)"
77+
]
78+
},
79+
{
80+
"cell_type": "code",
81+
"execution_count": null,
82+
"metadata": {},
83+
"outputs": [],
84+
"source": [
85+
"y_col.T @ y_left_null # column space is ortho to left null space"
86+
]
87+
},
88+
{
89+
"cell_type": "code",
90+
"execution_count": null,
91+
"metadata": {},
92+
"outputs": [],
93+
"source": [
94+
"# magnitudes of vectors\n",
95+
"np.sqrt(y_col.T @ y_col), np.sqrt(y_left_null.T @ y_left_null)"
96+
]
97+
},
98+
{
99+
"cell_type": "code",
100+
"execution_count": null,
101+
"metadata": {},
102+
"outputs": [],
103+
"source": [
104+
"X.T @ X # this is full rank -> invertible"
105+
]
106+
},
107+
{
108+
"cell_type": "code",
109+
"execution_count": null,
110+
"metadata": {},
111+
"outputs": [],
112+
"source": [
113+
"inv(X.T @ X)"
114+
]
115+
},
116+
{
117+
"cell_type": "code",
118+
"execution_count": null,
119+
"metadata": {},
120+
"outputs": [],
121+
"source": [
122+
"# left inverse for tall/thin, full column rank X\n",
123+
"Xli = inv(X.T @ X) @ X.T\n",
124+
"Xli"
125+
]
126+
},
127+
{
128+
"cell_type": "code",
129+
"execution_count": null,
130+
"metadata": {},
131+
"outputs": [],
132+
"source": [
133+
"# left inverse via SVD option 1 -> invert singular values & reverse space mapping: U -> V\n",
134+
"S = diagsvd(s, M, N)\n",
135+
"Sli = inv(S.T @ S) @ S.T\n",
136+
"Xli_svd_1 = V @ Sli @ U.T"
137+
]
138+
},
139+
{
140+
"cell_type": "code",
141+
"execution_count": null,
142+
"metadata": {},
143+
"outputs": [],
144+
"source": [
145+
"# left inverse via SVD option 2 -> invert singular values & reverse space mapping: U -> V\n",
146+
"# s / s^2 = 1 / s might be nicer seen here\n",
147+
"Xli_svd_2 = V @ diagsvd(s / s**2, N, M) @ U.T\n",
148+
"\n",
149+
"np.allclose(Xli_svd_2, Xli_svd_1), np.allclose(Xli, Xli_svd_1), np.allclose(Xli, pinv(X))"
150+
]
151+
},
152+
{
153+
"cell_type": "code",
154+
"execution_count": null,
155+
"metadata": {},
156+
"outputs": [],
157+
"source": [
158+
"theta_hat = Xli @ y # it is rarely called that way in this context, but: we actually train a model with this operation\n",
159+
"theta_hat # fitted / trained model parameters"
160+
]
161+
},
162+
{
163+
"cell_type": "code",
164+
"execution_count": null,
165+
"metadata": {},
166+
"outputs": [],
167+
"source": [
168+
"Xli @ y_col # we get same theta_hat if using only column space stuff of y "
169+
]
170+
},
171+
{
172+
"cell_type": "code",
173+
"execution_count": null,
174+
"metadata": {},
175+
"outputs": [],
176+
"source": [
177+
"Xli @ y_left_null # this must yield zero, as X cannot bring left null to row space"
178+
]
179+
},
180+
{
181+
"cell_type": "code",
182+
"execution_count": null,
183+
"metadata": {},
184+
"outputs": [],
185+
"source": [
186+
"y_hat = X @ theta_hat\n",
187+
"y_hat # == y_col"
188+
]
189+
},
190+
{
191+
"cell_type": "code",
192+
"execution_count": null,
193+
"metadata": {},
194+
"outputs": [],
195+
"source": [
196+
"e = y - y_hat # e == y_lns\n",
197+
"e, e.T @ e"
198+
]
199+
},
200+
{
201+
"cell_type": "code",
202+
"execution_count": null,
203+
"metadata": {},
204+
"outputs": [],
205+
"source": [
206+
"y_col.T @ e # column space is ortho to left null space"
207+
]
208+
},
209+
{
210+
"cell_type": "code",
211+
"execution_count": null,
212+
"metadata": {},
213+
"outputs": [],
214+
"source": [
215+
"# projection matrices\n",
216+
"\n",
217+
"P_col = X @ Xli\n",
218+
"P_col, P_col @ y, y_col"
219+
]
220+
},
221+
{
222+
"cell_type": "code",
223+
"execution_count": null,
224+
"metadata": {},
225+
"outputs": [],
226+
"source": [
227+
"# check projection in terms of SVD\n",
228+
"S @ Sli, np.allclose(U @ S @ Sli @ U.T, P_col)"
229+
]
230+
},
231+
{
232+
"cell_type": "code",
233+
"execution_count": null,
234+
"metadata": {},
235+
"outputs": [],
236+
"source": [
237+
"P_left_null = np.eye(M) - P_col\n",
238+
"P_left_null, P_left_null @ y, e"
239+
]
240+
},
241+
{
242+
"cell_type": "code",
243+
"execution_count": null,
244+
"metadata": {},
245+
"outputs": [],
246+
"source": [
247+
"P_row = Xli @ X # == always identity matrix for full column rank X\n",
248+
"P_row, P_row @ theta_hat"
249+
]
250+
},
251+
{
252+
"cell_type": "code",
253+
"execution_count": null,
254+
"metadata": {},
255+
"outputs": [],
256+
"source": [
257+
"# check projection in terms of SVD\n",
258+
"Sli @ S, np.allclose(V @ Sli @ S @ V.T, P_row)"
259+
]
260+
},
261+
{
262+
"cell_type": "code",
263+
"execution_count": null,
264+
"metadata": {},
265+
"outputs": [],
266+
"source": [
267+
"P_null = np.eye(N) - P_row # == always zero matrix for full column rank X\n",
268+
"P_null # null space is spanned only by zero vector"
269+
]
270+
},
271+
{
272+
"cell_type": "code",
273+
"execution_count": null,
274+
"metadata": {},
275+
"outputs": [],
276+
"source": [
277+
"plt.figure(figsize=(8,4))\n",
278+
"\n",
279+
"# residuals\n",
280+
"for m in range(M):\n",
281+
" plt.plot([X[m, 1], X[m, 1]],\n",
282+
" [y[m, 0], y_col[m, 0]], lw=3, label='error '+str(m+1))\n",
283+
"# data\n",
284+
"plt.plot(X[:,1], y, 'C4x',\n",
285+
" ms=10, mew=3,\n",
286+
" label='data')\n",
287+
"# fitted line\n",
288+
"plt.plot(X[:,1], theta_hat[0] * X[:,0] + theta_hat[1] * X[:,1], 'k', label='LS fit (interpolation)')\n",
289+
"x = np.linspace(0, 1, 10)\n",
290+
"plt.plot(x, theta_hat[0] + theta_hat[1] * x, 'C7:', label='LS fit (extrapolation)')\n",
291+
"x = np.linspace(4, 5, 10)\n",
292+
"plt.plot(x, theta_hat[0] + theta_hat[1] * x, 'C7:')\n",
293+
"\n",
294+
"plt.xticks(np.arange(6))\n",
295+
"plt.yticks(np.arange(11)-1)\n",
296+
"plt.xlim(0, 5)\n",
297+
"plt.ylim(-1, 9)\n",
298+
"plt.xlabel('feature x1')\n",
299+
"plt.ylabel('y')\n",
300+
"plt.title(r'min the sum of squared errors solves for $\\hat{\\theta}=[-1,2]^T$ -> intercept: -1, slope: +2')\n",
301+
"plt.legend()\n",
302+
"plt.grid(True)"
303+
]
304+
},
305+
{
306+
"cell_type": "markdown",
307+
"metadata": {},
308+
"source": [
309+
"## Copyright\n",
310+
"\n",
311+
"- the notebooks are provided as [Open Educational Resources](https://en.wikipedia.org/wiki/Open_educational_resources)\n",
312+
"- the text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/)\n",
313+
"- the code of the IPython examples is licensed under the [MIT license](https://opensource.org/licenses/MIT)\n",
314+
"- feel free to use the notebooks for your own purposes\n",
315+
"- please attribute the work as follows: *Frank Schultz, Data Driven Audio Signal Processing - A Tutorial Featuring Computational Examples, University of Rostock* ideally with relevant file(s), github URL https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise, commit number and/or version tag, year."
316+
]
317+
}
318+
],
319+
"metadata": {
320+
"kernelspec": {
321+
"display_name": "myddasp",
322+
"language": "python",
323+
"name": "python3"
324+
},
325+
"language_info": {
326+
"codemirror_mode": {
327+
"name": "ipython",
328+
"version": 3
329+
},
330+
"file_extension": ".py",
331+
"mimetype": "text/x-python",
332+
"name": "python",
333+
"nbconvert_exporter": "python",
334+
"pygments_lexer": "ipython3",
335+
"version": "3.12.3"
336+
}
337+
},
338+
"nbformat": 4,
339+
"nbformat_minor": 4
340+
}

0 commit comments

Comments
 (0)