\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}