163 lines
9.0 KiB
TeX
163 lines
9.0 KiB
TeX
|
\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}
|