Compare commits
No commits in common. "188ec4371b094060b20832676b876851e002f0b1" and "8f6a9e491dc01940ab2f56e02b684f180621d6ab" have entirely different histories.
188ec4371b
...
8f6a9e491d
Binary file not shown.
Binary file not shown.
@ -1,123 +0,0 @@
|
||||
\documentclass{article}
|
||||
\usepackage[utf8]{inputenc}
|
||||
\usepackage{biblatex}
|
||||
\addbibresource{library.bib}
|
||||
\usepackage{listings}
|
||||
\usepackage{amssymb}
|
||||
\usepackage{graphicx,amsmath}
|
||||
\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
|
||||
\usepackage{hyperref}
|
||||
\hypersetup{
|
||||
colorlinks=true,
|
||||
linkcolor=blue,
|
||||
filecolor=magenta,
|
||||
urlcolor=cyan,
|
||||
pdftitle={Overleaf Example},
|
||||
pdfpagemode=FullScreen,
|
||||
}
|
||||
|
||||
\title{Numerical Methods: Lecture 3. Projectors. Least squares problem. QR factorization.}
|
||||
\author{Konstantin Tikhonov}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\maketitle
|
||||
|
||||
\section{Suggested Reading}
|
||||
|
||||
\begin{itemize}
|
||||
\item Lectures 6-8, 10-11 of \cite{trefethen1997numerical}
|
||||
\item Lecture 8 of \cite{tyrtyshnikov2012brief}
|
||||
\end{itemize}
|
||||
|
||||
\section{Exercises}
|
||||
|
||||
Deadline: 4 Nov
|
||||
|
||||
\begin{enumerate}
|
||||
|
||||
\item (3) Consider the matrices:
|
||||
$$
|
||||
A=\quad\begin{bmatrix}
|
||||
1 & 0\\
|
||||
0 & 1\\
|
||||
1 & 0
|
||||
\end{bmatrix},\quad
|
||||
B=\quad\begin{bmatrix}
|
||||
1 & 2\\
|
||||
0 & 1\\
|
||||
1 & 0
|
||||
\end{bmatrix}.
|
||||
$$
|
||||
\begin{itemize}
|
||||
\item Derive orthogonal projectors on $\mathrm{range}(A)$ and $\mathrm{range}(B)$.
|
||||
\item Derive (on a piece of paper) QR decomposition of matrices $A$ and $B$.
|
||||
\end{itemize}
|
||||
|
||||
\item (5)
|
||||
Consider a particle of unit mass, which is prepared at $t=0$ at $x=0$ at rest $v=0$. The particle is exposed to piece--wise constant external force $f_i$ at $i-1< t \le i$, with $i=1,2,...,10$. Let $a=(x(t=10),v(t=10))$ be a vector composed of coordinate and velocity of a particle at $t=10$. Derive the matrix $A$ such that $a=Af$ (note that $A$ is of a shape $2\times 10$). Using (a numerical) SVD decomposition, evaluate $f$ of minimal norm such that $a=(1,0)$.
|
||||
|
||||
\item (5) Consider the function $f(x) = 10 \sin(x)$. Generate a dataset $D$ that will consist
|
||||
of $n = 7$ points drawn as follows. For each point randomly draw $x_i$ uniformly in $[0,6]$ and define $y_i = f(x_i) + \epsilon_i$,
|
||||
where $\epsilon_i$ are iid standard gaussian random numbers. Generate a sample dataset from this distribution, plot it together with the true function $f(x)$. Fit a linear $l(x) = w_0 + w_1x $ and a cubic $c(x) = w_0 + w_1x + w_2x^2 + w_3x^3$ models to D. Plot those models together with the dataset $D$.
|
||||
|
||||
\item (7) Download the \href{https://www.dropbox.com/s/qgz1x67t10fd7hf/data.npz?dl=0}{file} with matrices $A$ and $C$ (an image and a filter). Open it as follows: \lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
\lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
with np.load('data.npz') as data:
|
||||
A, C = data['A'], data['C']
|
||||
\end{lstlisting}
|
||||
It is convenient to order the matrix $A$ into a column vector $a$:
|
||||
\lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
\lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
def mat2vec(A):
|
||||
A = np.flipud(A)
|
||||
a = np.reshape(A, np.prod(A.shape))
|
||||
return a
|
||||
\end{lstlisting}
|
||||
with inverse transform, from vector $a$ to matrix $A$ given by
|
||||
\lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
\lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
def vec2mat(a, shape):
|
||||
A = np.reshape(a, shape)
|
||||
A = np.flipud(A)
|
||||
return A
|
||||
\end{lstlisting}
|
||||
The image, stored in the matrix $A$ is obtained from certain original image $A_0$ via convoluting it with the filter $C$ and adding some noise. The filter $C$ blurs an image, simultaneously increasing its size from $16\times 51$ to $25\times 60$. With the use of associated vectors $a$ and $a_0$, one may write
|
||||
$$
|
||||
a_0\to a = C a_0 + \epsilon,
|
||||
$$
|
||||
where $\epsilon$ is a vector of iid Gaussian random numbers. Your task will be to recover an original image $A_0$, being supplied by the image $A$ and the filter $C$.
|
||||
\begin{itemize}
|
||||
\item Plot the image $A$.
|
||||
\item Explore how the filter $C$ acts on images.
|
||||
\item A naive way to recover $A_0$ from $A$ would be to solve $a = C a_0$ for $a_0$. Is this system under-- or over--determined? Using SVD of the filter matrix $C$, evaluate $a_0$ and plot the corresponding $A_0$.
|
||||
\item In order to improve the result, consider keeping certain fraction of singular values of $C$. Choose a value delivering the best recovery quality.
|
||||
\end{itemize}
|
||||
\item (7) Consider the problem
|
||||
$$
|
||||
\mathrm{minimize\quad}\Vert Ax-b \Vert_2 \mathrm{\quad subject\;to\quad}Cx=0\mathrm{\quad with\;respect\;to\quad}x.
|
||||
$$
|
||||
Using the method of Lagrange multipliers, and assuming $A^TA$ to be invertible, derive explicit expression for optimal $x$.
|
||||
|
||||
\item (20) Here we will consider the problem of localization of points in a 2D plane. Consider $n$ points, for which we have the \emph{approximate} locations $r_i=\left(x_i, y_i\right)$. We measure $k$ angles between certain points: $\theta_{ijk}=\angle(r_k-r_i, r_j-r_i)$. Our goal is to use the results of the measurements to improve the estimation of the locations $r_i$.
|
||||
|
||||
To be specific, consider $n=3$ points, for which we have approximate locations $r_1=(-1,0),\;r_2=(0, 1),\;r_3=(1,0)$ and $k=1$ measurement $\theta_{123}=9\pi/40$. Clearly, our estimates $r_{1,2,3}$ are not consistent with the measured angle and we have to adjust the estimate: $r_i\to \bar{r}_i= r_i+dr_i$ where $dr_i$ should be found from the condition
|
||||
$$
|
||||
(\bar{r}_3-\bar{r}_1)\cdot(\bar{r}_3-\bar{r}_1) = |\bar{r}_3-\bar{r}_1||\bar{r}_3-\bar{r}_1|\cos\theta_{123},
|
||||
$$
|
||||
which can be linearized in $dr_i$, assuming this correction will end up small. In this approach, one constructs a single (in general, $k$) equation for six (in general, $2n$) variables, so the system will typically by overdetermined. We can consider this system in the least squares sense, which amounts to determining the smallest correction to all $r_i$ which makes the updated locations consistent with observations. In the particular numerical example above, one may find $dr_1=(-h, 0),\;dr_2=(h, -h),\; dr_3=(0, h)$ where $h=\pi/80\approx 0.04$.
|
||||
|
||||
Your task is to write the code, which will accept the current estimate of the positions $r_i$ ($n\times 2$, float) and measurement results $\theta_{ijk}$, which are specified by i) indices of points ($k\times 3$, int) and ii) angles ($k$, float); and will output the derived correction to the point positions $dr_i$ ($n\times 2$, float). You can use the numerical example in this exercise to test your code.
|
||||
|
||||
\end{enumerate}
|
||||
\printbibliography
|
||||
\end{document}
|
||||
@ -1,162 +0,0 @@
|
||||
\documentclass[prb, notitlepage, aps, 11pt]{revtex4-2}%
|
||||
\usepackage[utf8]{inputenc}
|
||||
\usepackage[T2A]{fontenc}
|
||||
\usepackage[english, russian]{babel}
|
||||
%\usepackage{listings}
|
||||
\usepackage{amssymb}
|
||||
\usepackage{graphicx,amsmath}
|
||||
\usepackage{enumitem}
|
||||
\usepackage{nicefrac}
|
||||
\usepackage{amsmath}
|
||||
\usepackage{graphicx}
|
||||
\usepackage{amsfonts}
|
||||
\usepackage{comment}
|
||||
\usepackage{bm}
|
||||
\usepackage{delimset}
|
||||
\usepackage[pdftitle = a]{hyperref}
|
||||
\usepackage{datetime}
|
||||
\usepackage{minted}
|
||||
\usemintedstyle{friendly}
|
||||
\usepackage[a]{esvect}
|
||||
\hypersetup{
|
||||
colorlinks=true,
|
||||
linkcolor=blue,
|
||||
filecolor=magenta,
|
||||
urlcolor=cyan,
|
||||
pdfpagemode=FullScreen,
|
||||
}
|
||||
\usepackage{microtype}
|
||||
|
||||
\newcommand{\framesep}{0.6em}
|
||||
\BeforeBeginEnvironment{minted}{\vspace{-1.6em}}
|
||||
\AfterEndEnvironment{minted}{\vspace{-0.5em}}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\begin{center}
|
||||
версия от \today \quad \currenttime
|
||||
\end{center}
|
||||
|
||||
\title{\texorpdfstring{
|
||||
Численные методы, осень 2022\\
|
||||
Задание 3 [Проекторы. Задача наименьших квадратов. QR-разложение.] \\
|
||||
Всего баллов: 47 \ Срок сдачи: 4 ноября
|
||||
}{}
|
||||
}
|
||||
|
||||
\maketitle
|
||||
|
||||
\section*{Рекомендованная литература}
|
||||
|
||||
\begin{itemize}
|
||||
\item Лекции 6--8, 10--11 из \cite{trefethen1997numerical}
|
||||
\item Лекция 8 из \cite{tyrtyshnikov2012brief}
|
||||
\end{itemize}
|
||||
|
||||
\section*{Упражнения}
|
||||
|
||||
\begin{enumerate}
|
||||
|
||||
\item (3) Даны матрицы:
|
||||
$$
|
||||
A=\quad\begin{bmatrix}
|
||||
1 & 0\\
|
||||
0 & 1\\
|
||||
1 & 0
|
||||
\end{bmatrix},\quad
|
||||
B=\quad\begin{bmatrix}
|
||||
1 & 2\\
|
||||
0 & 1\\
|
||||
1 & 0
|
||||
\end{bmatrix}
|
||||
$$
|
||||
\begin{itemize}
|
||||
\item Получите ортогональные проекторы на пространства $\operatorname{colsp}(A)$ и $\operatorname{colsp}(B)$.
|
||||
\item Вычислите (вручную) QR-разложения матриц $A$ and $B$.
|
||||
\end{itemize}
|
||||
|
||||
\item (5)
|
||||
Частица единичной массы в начальный момент времени покоится в точке $x = 0$. Затем она подвергается воздействию кусочно-постоянной внешней силы $f_i \;\; (i = 1 \ldots 10)$.
|
||||
Пусть $a = \eval{(x,\: v)}_{t=10}$ --- вектор, состоящий из координаты и скорости частицы в конечный момент времени. Найдите матрицу $A$ такую, что $a = Af$. Заметьте, что она будет иметь размер $2\times 10$. С помощью SVD-разложения численно найдите вектор $f$ наименьшей нормы, при котором $a = (1,\, 0)$.
|
||||
|
||||
\item (5)
|
||||
Создайте датасет следующим образом: выберите $n = 7$ точек $x_i$ из равномерного распределения на отрезке $[0,\,6]$. Затем вычислите $y_i = f(x_i) + \epsilon_i$, где $f(x) = 10 \sin(x)$, а $\epsilon_i$~---~независимые стандартные нормальные случайные величины. Изобразите на одном графике точки датасета и функцию $f$.
|
||||
|
||||
Аппроксимируйте датасет с помощью линейной $l(x) = w_0 + w_1 x$ и кубической $c(x) = w_0 + w_1x + w_2x^2 + w_3x^3$ функций. Изобразите на одном графике получившиеся функции и исходный датасет.
|
||||
|
||||
%\pagebreak
|
||||
|
||||
\item (7)
|
||||
Загрузите \href{https://www.dropbox.com/s/qgz1x67t10fd7hf/data.npz?dl=0}{файл} с матрицами $A$ и $C$ (изображение и фильтр). Откройте его:
|
||||
|
||||
\begin{minted}[frame = lines, framesep = \framesep]{python}
|
||||
with np.load('data.npz') as data:
|
||||
A, C = data['A'], data['C']
|
||||
\end{minted}
|
||||
|
||||
Матрицу $A$ полезно преобразовать в вектор $a$:
|
||||
\begin{minted}[frame = lines, framesep = \framesep]{python}
|
||||
def mat2vec(A):
|
||||
A = np.flipud(A)
|
||||
a = np.reshape(A, np.prod(A.shape))
|
||||
return a
|
||||
\end{minted}
|
||||
|
||||
Обратное преобразование из вектора $a$ в матрицу $A$:
|
||||
\begin{minted}[frame = lines, framesep = \framesep]{python}
|
||||
def vec2mat(a, shape):
|
||||
A = np.reshape(a, shape)
|
||||
A = np.flipud(A)
|
||||
return A
|
||||
\end{minted}
|
||||
|
||||
Изображение, хранящееся в матрице $A$, было получено из исходного изображения $A_0$ с помощью свёртки с фильтром $C$ и добавления шума. Записывая матрицы как векторы описанным выше способом, мы можем написать
|
||||
$$
|
||||
a_0\to a = C a_0 + \epsilon,
|
||||
$$
|
||||
где $\epsilon$ --- вектор из независимых одинаково распределённых гауссовских случайных величин.
|
||||
Фильтр $C$ размывает изображение, увеличивая при этом его размер с $16 \times 51$ до $25 \times 60$.
|
||||
Вашей задачей будет восстановить исходное изображение $A_0$ по размытому изображению $A$ и фильтру $C$.
|
||||
%
|
||||
\begin{itemize}
|
||||
\item Постройте изображение $A$.
|
||||
|
||||
\item Исследуйте, как фильтр $C$ меняет изображения.
|
||||
|
||||
\item Наивный способ восстановить $A_0$ --- это найти $a_0$ из системы $a = Ca_0$. Является ли эта система недо- или переопределённой? Вычислите $a_0$ с помощью сингулярного разложения матрицы $C$ и постройте получившееся изображение.
|
||||
|
||||
\item Чтобы улучшить результат, попробуйте использовать только часть сингулярных чисел матрицы $C$. Выберите наиболее удачное количество.
|
||||
\end{itemize}
|
||||
%
|
||||
\item (7) Рассмотрим задачу оптимизации
|
||||
$$
|
||||
\min_x\, \norm!{Ax - b} \quad\text{при условии}\quad Cx = 0
|
||||
$$
|
||||
Предполагая, что матрица $A^T A$ обратима, воспользуйтесь методом множителей Лагранжа и получите выражение для точки минимума $x$.
|
||||
|
||||
\item (20) Эта задача посвящена локализации точек на плоскости. Рассмотрим $n$ точек с \emph{неточными} координатами: $A_i = \brk{x_i,\, y_i}$. Измерим $m$ углов между некоторыми из них: $\theta_{ijk} = \angle\, A_i A_j A_k$. Наша цель --- используя результаты измерений, уточнить координаты точек.
|
||||
|
||||
В качестве примера рассмотрим $n = 3$ точки с приближенными координатами: $A_1 = (1,\, 0),\ A_2 = (-1,\, 0),\ A_3 = (0,\, 1)$ --- и один измеренный угол $\theta_{123} = 9\pi/40$. Так как наши оценки координат не согласуются с измеренным углом, мы должны их скорректировать:
|
||||
$$
|
||||
x_i \to x_i + dx_i, \quad y_i \to y_i + dy_i,
|
||||
$$
|
||||
где $dx_i$ и $dy_i$ находятся из условия
|
||||
$$
|
||||
\vv{A_1 A}_2 \cdot \vv{A_3 A}_2 =
|
||||
\abs{A_1 A_2} \cdot \abs{A_3 A_2} \cdot \cos(\theta_{123})
|
||||
$$
|
||||
Это условие может быть линеаризовано в предположении, что поправки окажутся малыми.
|
||||
При таком подходе мы конструируем $m = 1$ уравнений для $2n = 6$ переменных. В общем случае число измеренных углов $m$ может быть большим, и система окажется переопределённой. Но мы можем рассматривать её в смысле метода наименьших квадратов, т.е. пытаться определить наименьшие поправки к координатам, которые делали бы их согласованными с измерениями. В рассмотренном примере можно найти решение
|
||||
$$
|
||||
dr_1 = (dx_1,\, dy_1) = (0, h), \quad
|
||||
dr_2=(-h,\, 0), \quad
|
||||
dr_3=(h,\, -h), \qquad
|
||||
\text{ где } h = \pi/80 \approx 0.04
|
||||
$$
|
||||
Ваша задача --- написать код, который принимает текущие оценки координат точек $(n \times 2)$ и измеренные углы $\theta_{ijk}$ ($m$ значений и $m \times 3$ индексов). Возвращать он должен найденные поправки $(n \times 2)$.
|
||||
Вы можете использовать пример из этого упражнения, чтобы протестировать свой код.
|
||||
|
||||
\end{enumerate}
|
||||
|
||||
\bibliography{library.bib}
|
||||
\end{document}
|
||||
@ -1,118 +0,0 @@
|
||||
@book{trefethen1997numerical,
|
||||
title={Numerical linear algebra},
|
||||
author={Trefethen, Lloyd N and Bau III, David},
|
||||
volume={50},
|
||||
year={1997},
|
||||
publisher={Siam}
|
||||
}
|
||||
@book{robert2013monte,
|
||||
title={Monte Carlo statistical methods},
|
||||
author={Robert, Christian and Casella, George},
|
||||
year={2013},
|
||||
publisher={Springer Science \& Business Media}
|
||||
}
|
||||
@book{murphy2012machine,
|
||||
title={Machine learning: a probabilistic perspective},
|
||||
author={Murphy, Kevin P},
|
||||
year={2012},
|
||||
publisher={MIT press}
|
||||
}
|
||||
@book{boyd2004convex,
|
||||
title={Convex optimization},
|
||||
author={Boyd, Stephen and Boyd, Stephen P and Vandenberghe, Lieven},
|
||||
year={2004},
|
||||
publisher={Cambridge university press}
|
||||
}
|
||||
@book{trefethen2019approximation,
|
||||
title={Approximation Theory and Approximation Practice, Extended Edition},
|
||||
author={Trefethen, Lloyd N},
|
||||
year={2019},
|
||||
publisher={SIAM}
|
||||
}
|
||||
@book{devroye:1986,
|
||||
author = {Devroye, Luc},
|
||||
date = {1986)},
|
||||
description = {Non-Uniform Random Variate Generation},
|
||||
keywords = {algorithms generation random simulation},
|
||||
location = {New York},
|
||||
publisher = {Springer-Verlag},
|
||||
title = {Non-Uniform Random Variate Generation(originally published with},
|
||||
year = 1986
|
||||
}
|
||||
@book{williams2006gaussian,
|
||||
title={Gaussian processes for machine learning},
|
||||
author={Williams, Christopher K and Rasmussen, Carl Edward},
|
||||
volume={2},
|
||||
number={3},
|
||||
year={2006},
|
||||
publisher={MIT press Cambridge, MA}
|
||||
}
|
||||
@book{hairer1993solving,
|
||||
title={Solving ordinary differential equations. 1, Nonstiff problems},
|
||||
author={Hairer, Ernst and N{\o}rsett, Syvert P and Wanner, Gerhard},
|
||||
year={1993},
|
||||
publisher={Springer-Vlg}
|
||||
}
|
||||
@book{hairer1993solving2,
|
||||
title={Solving ordinary differential equations. 2, Stiff and differential-algebraic problems},
|
||||
author={Hairer, Ernst and N{\o}rsett, Syvert P and Wanner, Gerhard},
|
||||
year={1993},
|
||||
publisher={Springer-Vlg}
|
||||
}
|
||||
@book{tyrtyshnikov2012brief,
|
||||
title={A brief introduction to numerical analysis},
|
||||
author={Tyrtyshnikov, Eugene E},
|
||||
year={2012},
|
||||
publisher={Springer Science \& Business Media}
|
||||
}
|
||||
@book{amosov2003,
|
||||
title={Numerical Methods for Engineers},
|
||||
author={Amosov, AA and Dubinskiy YuA and Kopchenova, NV},
|
||||
year={2003},
|
||||
publisher={Izdatelstvo MEI}
|
||||
}
|
||||
@article{arbenz2012lecture,
|
||||
title={Lecture notes on solving large scale eigenvalue problems},
|
||||
author={Arbenz, Peter and Kressner, Daniel and Z{\"u}rich, DME},
|
||||
journal={D-MATH, EHT Zurich},
|
||||
volume={2},
|
||||
year={2012}
|
||||
}
|
||||
@article{trefethen1996finite,
|
||||
title={Finite difference and spectral methods for ordinary and partial differential equations},
|
||||
author={Trefethen, Lloyd Nicholas},
|
||||
year={1996},
|
||||
publisher={Cornell University-Department of Computer Science and Center for Applied~…}
|
||||
}
|
||||
@book{boyd2018introduction,
|
||||
title={Introduction to applied linear algebra: vectors, matrices, and least squares},
|
||||
author={Boyd, Stephen and Vandenberghe, Lieven},
|
||||
year={2018},
|
||||
publisher={Cambridge university press}
|
||||
}
|
||||
@article{halko2011finding,
|
||||
title={Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions},
|
||||
author={Halko, Nathan and Martinsson, Per-Gunnar and Tropp, Joel A},
|
||||
journal={SIAM review},
|
||||
volume={53},
|
||||
number={2},
|
||||
pages={217--288},
|
||||
year={2011},
|
||||
publisher={SIAM}
|
||||
}
|
||||
@book{demmel1997applied,
|
||||
title={Applied numerical linear algebra},
|
||||
author={Demmel, James W},
|
||||
year={1997},
|
||||
publisher={SIAM}
|
||||
}
|
||||
@article{dahlquist198533,
|
||||
title={33 years of numerical instability, Part I},
|
||||
author={Dahlquist, Germund},
|
||||
journal={BIT Numerical Mathematics},
|
||||
volume={25},
|
||||
number={1},
|
||||
pages={188--204},
|
||||
year={1985},
|
||||
publisher={Springer}
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user