\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \ifxetex
    \usepackage{mathspec}
    \makeatletter % undo the wrong changes made by mathspec
    \let\RequirePackage\original@RequirePackage
    \let\usepackage\RequirePackage
    \makeatother
  \else
    \usepackage{fontspec}
  \fi
  \defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage{hyperref}
\hypersetup{unicode=true,
            pdftitle={Homework 5},
            pdfborder={0 0 0},
            breaklinks=true}
\urlstyle{same}  % don't use monospace font for urls
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}
% Redefines (sub)paragraphs to behave more like sections
\ifx\paragraph\undefined\else
\let\oldparagraph\paragraph
\renewcommand{\paragraph}[1]{\oldparagraph{#1}\mbox{}}
\fi
\ifx\subparagraph\undefined\else
\let\oldsubparagraph\subparagraph
\renewcommand{\subparagraph}[1]{\oldsubparagraph{#1}\mbox{}}
\fi

\title{Homework 5}

\input{preamble.tex}
\title{Homework}



\begin{document}
\maketitle

\hypertarget{homework-5}{%
\section{Homework 5}\label{homework-5}}

Please answer the following questions in complete sentences in a clearly
prepared manuscript and submit the solution by the due date on
Gradescope, around 4am on Monday November 20nd

Remember that this is a graduate class. There may be elements of the
problem statements that require you to fill in appropriate assumptions.
You are also responsible for determining what evidence to include. An
answer alone is rarely sufficient, but neither is an overly verbose
description required. Use your judgement to focus your discussion on the
most interesting pieces. The answer to ``should I include `something' in
my solution?'' will almost always be: Yes, if you think it helps support
your answer.

\hypertarget{problem-0-homework-checklist}{%
\subsection{Problem 0: Homework
checklist}\label{problem-0-homework-checklist}}

\begin{itemize}
\item
  Please identify anyone, whether or not they are in the class, with
  whom you discussed your homework. This problem is worth 1 point, but
  on a multiplicative scale.
\item
  Make sure you have included your source-code and prepared your
  solution according to the most recent Piazza note on homework
  submissions.
\end{itemize}

\hypertarget{problem-1-updating-factorizations-of-linear-systems-of-equations}{%
\subsection{Problem 1: Updating factorizations of linear systems of
equations}\label{problem-1-updating-factorizations-of-linear-systems-of-equations}}

In class, we showed how to solve \((\mA + \vu \vv^T) \vx = \vb\) when
given a fast factorization method to solve \(\mA \vx = \vb\). In this
problem, we will address the question of how to update the factorization
itself. Suppose that we are given a Cholesky factorization of
\(\mA = \mL \mD \mL^T\) as we saw in class. Show how to update this
factorization to produce \(\mL' \mD' \mL'^T = (\mA + \vu \vu^T)\). Your
algorithm should do no more than \(O(n^2)\) work.

\hypertarget{problem-2-matrix-functions-in-practice}{%
\subsection{Problem 2: Matrix functions in
practice}\label{problem-2-matrix-functions-in-practice}}

In general, we don't compute Matrix functions by taking their
eigenvalues and applying a function to their eigenvalues -- this is just
a nice theory that works to explain what's going on in the simple case.

Instead, we create customized algorithms to compute a function of a
matrix. One such algorithm that computes the cosine of a matrix was
proposed by Serbin and Blalack and analyzed by Higham and Smith in 2003
\url{https://login.ezproxy.lib.purdue.edu/login?url=https://link.springer.com/article/10.1023/A:1026152731904}

The idea is that use the cosine double angle formula:
\[ \cos(2\mA) = 2 \cos(\mA)^2 - \mI \] combined with some way of
estimating \(\mC_0 = \cos(2^{-m} \mA)\). Then, we can compute
\(\cos(2^{-m+1} \mA) = 2 \mC_0^2 - \mI\). Also, for \(\cos(2^{-m} \mA)\)
then we should be able to use a simple formula to get a good enough
approximation.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Show that the double angle formula holds for symmetric matrices
  \(\mA\).
\item
  In the paper, they show that if \(\normof{\mA}_\infty \le 1\), then
  there is a specific type of approximation of \(\cos(\mA)\) called a
  Pad\textquotesingle\{e\} approxiomation that gives an answer down to
  machine precision. The specific approximation is
  \[ \cos(x) \approx r_{88}(x) = 
  \frac{1 - 
    \frac{260735}{545628} x^2 + 
    \frac{4375409}{141863280} x^4 -
    \frac{7696415}{13108167072} x^6 + 
    \frac{80737373}{23594700729600} x^8 }{ 
   1 + \frac{12079}{545628} x^2 + 
     \frac{34709}{141863280} x^4 + 
     \frac{109247}{65540835360} x^6 + 
     \frac{11321}{1814976979200} x^8
  }.  \]

  When applied to a matrix with \(\normof{\mA}_{\infty} \le 1\), this
  approximation guarantees
  \[ \frac{\normof{\cos(\mA) - r_{88}(\mA)}_{\infty}}{\normof{\cos(\mA)}_{\infty}} \le 3.26 \times 10^{-16}. \]

  Note that this is a rational function \(r_{88}(x) = p(x) / q(x)\) for
  some polynomials \(p(x)\) and \(q(x)\). One way to evaluate this for a
  matrix is compute \(p(\mA) \cdot q(\mA)^{-1}\).

  Based on \(r_{88}\), let's create values \(\pi\) and \(\mu\) so that
  \[ p(x) = \pi_0 + \pi_2 x^2 + \pi_4 x^4 + \pi_6 x^6 + \pi_8 x^8 \]
  \[ q(x) = \mu_0 + \mu_2 x^2 + \mu_4 x^4 + \mu_6 x^6 + \mu_8 x^8 \]
  Then to evaluate \(r_{88}(\mA)\) we have the algorithm.

\begin{verbatim}
     A2 = A*A
     A4 = A2*A2 
     A6 = A4 * A2
     A8 = A4 * A4
     P = pi0*I + pi2*A2 + pi4*A4 + pi6*A6 + pi8*A8
     Q = mu0*I + mu2*A2 + mu4*A4 + mu6*A6 + mu8*A8
     R = Q\P
\end{verbatim}

  Implement this algorithm and show that for matrices where
  \(\normof{\mA}_\infty \le 1\) we have the given bound on accuracy.
  (\textbf{Instruction note: While I have tried to be careful to check
  the accuracy of the coefficients, we will verify this ourselves and
  let you know if there are issues and note any correction to the
  problem.})

  You may assume that calling \(\cos(\mA)\) in Julia is accurate.
\item
  Show how to determine the smallest \(m\) such that
  \(\normof{2^{-m} \mA}_\infty \le 1\).
\item
  Implement the Serbin/Blalack/Higham/Smith procedure to compute
  \(\cos(\mA)\), i.e.~

\begin{verbatim}
 compute m such that ||2^{-m} A||_oo <= 1 
 compute R = r_88(2^{-m} A)
 for i=1:m 
     R = 2*R*R - I 
 end
 return R 
\end{verbatim}

  For a 1000 by 1000 matrix, compare the runtime to compute
  \(\cos(\mA)\) via a spectral eigenvalue decomposition.
\end{enumerate}

\hypertarget{problem-3-weighted-orthogonality-for-svd.}{%
\subsection{Problem 3: Weighted orthogonality for
SVD.}\label{problem-3-weighted-orthogonality-for-svd.}}

One possible weighted generalization of the SVD involves producing a
factorization \[ \mA = \mU \mSigma \mV^T \] where \(\mU^T \mR \mU=\mI\)
and \(\mV^T \mS \mV = \mI\), where \$\mR \$ and \(\mS\) are symmetric
positive definite matrices.

A common theme in solving more advanced problems is converting or
translating them into problems that we know how to solve.

Show how to use a standard SVD computation to produce this weighted SVD
factorization.

\hypertarget{problem-4-implementing-a-rank-k-solution-update.}{%
\subsection{Problem 4: Implementing a rank-k solution
update.}\label{problem-4-implementing-a-rank-k-solution-update.}}

Given \(\vx\) and an LU factorization of \(\mA\), then updating the
solution of a linear system \(\mA \vx = \vb\) to a new solution
\((\mA + \mU \mV^T) \vy = \vb\) can be done fairly efficiently.

Write code or fairly detailed pseudo-code to do this. This uses the
Julia factorization \texttt{F} which enables you to solve systems with
\(\mA\) without recomputing the factorization.

\begin{verbatim}
function update(x::Vector, b::Vector, F::LUFact, U::Matrix, V::Matrix)
end
\end{verbatim}

\hypertarget{problem-5-adding-and-deleting-an-equation}{%
\subsection{Problem 5: Adding and deleting an
equation}\label{problem-5-adding-and-deleting-an-equation}}

Recall that a linear system represents the simultaneously solution of a
set of linear equations \[ \va_1^T \vx = b_1 \] \[ \va_2^T \vx = b_2 \]
\[ \vdots \] \[ \va_n^T \vx = b_n. \] Consider a set of \(n\) equations
and \(n\) unknowns with a unique solution for any possible set of values
\(b_1, \ldots, b_n\).

Let \(\texttt{alg}(\vb)\) be an algorithm to solve for \(\vx\) when
given \(b_1, \ldots b_n\).

Show how to use \(\texttt{alg}\) in order to solve for \(\vy\) in the
system of equations \[ \vc_1^T \vy = d_1 \] \[ \va_2^T \vy = b_2 \]
\[ \vdots \] \[ \va_n^T \vy = b_n. \] That is, we deleted the equation
with \(\va_1\) and added a new equation with \(\vc_1\) instead. If you
use \texttt{alg} too many times, you may lose points. Be efficient!

(Hint: this is essentially a special case of other problems on this
homework or other things we've seen in class.)

\hypertarget{problem-6-the-krylov-subspace-and-eigenvalue-algorithms}{%
\subsection{Problem 6: The Krylov subspace and eigenvalue
algorithms}\label{problem-6-the-krylov-subspace-and-eigenvalue-algorithms}}

In class, we showed that the Krylov subspace:
\[ \mathbb{K}_k(\mA, \vb) = \text{span}(\vb, \mA \vb, \ldots, \mA^{k} \vb) \]
arises when solving \(\mA \vx = \vb\) via the simple Neumann series
method. Here, we note that it also arises in the power-method to
identify the largest eigenvalue of a matrix. Note that the simple
monomial basis for the Krylov subspace is:
\[ \mX_k = \bmat{ \vb & \mA \vb & \ldots & \mA^{k} \vb} \] In this case,
the power-method uses the vector \(\mX_k \ve_k\) as the estimate of the
largest eigenvalue.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  State an algorithm to search the Krylov subspace for the largest
  eigenvalue of a symmetric matrix using the fact that the largest
  eigenvalue solves the problem:
  \[ \lambda_{\max} = \MAXone{}{\vx^T \mA \vx}{\normof{\vx} = 1.}  \]
  (Note, that you are welcome to use ideas that we discuss in subsequent
  lectures, but you don't have to. Everything can be done with the
  material presented through lecture 21.) Note that you may use any
  routine you want on a \((k+1)\)-by-\((k+1)\) matrix.
\item
  Implement your algorithm and compare the eigenvalue estimates for a
  few small matrices. (Hint, a good idea here would be to show
  convergence plots for the error in the largest eigenvalue compared
  with iterations -- usually on a log-y-scale.)
\item
  Comment on any issues that arise in your implementation.
\end{enumerate}

\end{document}
