numerics-2022/Task 3/tex/Assignment_3_ru.tex

171 lines
9.6 KiB
TeX
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

\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 \ Срок сдачи (после переноса): 11 ноября
}{}
}
\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}_2 \quad\text{при условии}\quad Cx = \vv 0
$$
Предполагая, что матрица $A^T A$ обратима, воспользуйтесь методом множителей Лагранжа и получите выражение для точки минимума $x$.
\item (20) Эта задача посвящена локализации точек на плоскости. Рассмотрим $n$ точек с \emph{неточными} координатами: $A_i = \brk{x_i,\, y_i}$. Измерим $m$ углов между некоторыми из них: $\theta_{ijk} = \angle \brk!{\vv{A_i A}_j,\ \vv{A_i A}_k}$. Наша цель --- используя результаты измерений, уточнить координаты точек.
В качестве примера рассмотрим $n = 3$ точки с приближенными координатами: $A_1 = (-1,\, 0),\ A_2 = (0,\, 1),\ A_3 = (1,\, 0)$ --- и один измеренный угол $\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_1 A}_3 =
\abs{A_1 A_2} \cdot \abs{A_1 A_3} \cdot \cos(\theta_{123})
$$
Это условие может быть линеаризовано в предположении, что поправки окажутся малыми.
При таком подходе мы конструируем $m = 1$ уравнений для $2n = 6$ переменных. Система обычно оказывается недоопределённой. Но мы можем рассматривать её в смысле метода наименьших квадратов, т.е. пытаться определить наименьшие поправки к координатам, которые делали бы их согласованными с измерениями. В рассмотренном примере можно найти решение
$$
dr_1 = (dx_1,\, dy_1) = (-h,\, 0), \quad
dr_2=(h,\, -h), \quad
dr_3=(0,\, h), \qquad
\text{ где } h \approx \pi/80 \approx 0.04
$$
Ваша задача --- написать код, который принимает текущие оценки координат точек $(n \times 2)$ и измеренные углы $\theta_{ijk}$ ($m$ значений и $m \times 3$ индексов). Возвращать он должен найденные поправки $(n \times 2)$.
Чтобы протестировать свой код, вы можете использовать \href{https://disk.yandex.ru/d/kNjRT_r3S2xk2A}{пример} из этого упражнения, а также другой \href{https://disk.yandex.ru/d/7aGeQvikiUt1AA}{пример} с б\'{о}льшим числом точек. Загрузить данные можно следующим образом:
\begin{minted}[frame = lines, framesep = \framesep]{python}
with np.load('./data_1.npz') as data:
r = data['r'] # исходные точки
p = data['p'] - 1 # индексы измеренных углов
theta = data['theta'] # углы в радианах
dr = data['dr'] # предполагаемый ответ
\end{minted}
\end{enumerate}
\bibliography{library.bib}
\end{document}