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

171 lines
9.6 KiB
TeX
Raw Permalink Normal View History

2022-10-15 00:00:42 +03:00
\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-разложение.] \\
2022-11-02 16:10:21 +03:00
Всего баллов: 47 \ Срок сдачи (после переноса): 11 ноября
2022-10-15 00:00:42 +03:00
}{}
}
\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)$.
2022-10-19 16:41:34 +03:00
Пусть $a = \eval{(x,\: v)}_{t=10}$ --- вектор, состоящий из координаты и скорости частицы в конечный момент времени. Найдите матрицу $A$ такую, что $a = Af$. Заметьте, что она будет иметь размер $2\times 10$. С помощью SVD-разложения численно найдите вектор $f$ наименьшей евклидовой нормы, при котором $a = (1,\, 0)$.
2022-10-15 00:00:42 +03:00
\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) Рассмотрим задачу оптимизации
$$
2022-10-19 16:41:34 +03:00
\min_x\, \norm!{Ax - b}_2 \quad\text{при условии}\quad Cx = \vv 0
2022-10-15 00:00:42 +03:00
$$
Предполагая, что матрица $A^T A$ обратима, воспользуйтесь методом множителей Лагранжа и получите выражение для точки минимума $x$.
2022-11-02 16:10:21 +03:00
\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}$. Наша цель --- используя результаты измерений, уточнить координаты точек.
2022-10-15 00:00:42 +03:00
2022-11-02 16:10:21 +03:00
В качестве примера рассмотрим $n = 3$ точки с приближенными координатами: $A_1 = (-1,\, 0),\ A_2 = (0,\, 1),\ A_3 = (1,\, 0)$ --- и один измеренный угол $\theta_{123} = 9\pi/40$. Так как наши оценки координат не согласуются с измеренным углом, мы должны их скорректировать:
2022-10-15 00:00:42 +03:00
$$
x_i \to x_i + dx_i, \quad y_i \to y_i + dy_i,
$$
где $dx_i$ и $dy_i$ находятся из условия
$$
2022-11-02 16:10:21 +03:00
\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})
2022-10-15 00:00:42 +03:00
$$
Это условие может быть линеаризовано в предположении, что поправки окажутся малыми.
2022-10-19 16:41:34 +03:00
При таком подходе мы конструируем $m = 1$ уравнений для $2n = 6$ переменных. Система обычно оказывается недоопределённой. Но мы можем рассматривать её в смысле метода наименьших квадратов, т.е. пытаться определить наименьшие поправки к координатам, которые делали бы их согласованными с измерениями. В рассмотренном примере можно найти решение
2022-10-15 00:00:42 +03:00
$$
2022-11-02 16:10:21 +03:00
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
2022-10-15 00:00:42 +03:00
$$
Ваша задача --- написать код, который принимает текущие оценки координат точек $(n \times 2)$ и измеренные углы $\theta_{ijk}$ ($m$ значений и $m \times 3$ индексов). Возвращать он должен найденные поправки $(n \times 2)$.
2022-11-02 16:10:21 +03:00
Чтобы протестировать свой код, вы можете использовать \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}
2022-10-15 00:00:42 +03:00
\end{enumerate}
\bibliography{library.bib}
\end{document}