\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 3b},
            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 3b}

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



\begin{document}
\maketitle

\section{Homework 3b}\label{homework-3b}

Please answer the following questions in complete sentences in a clearly
prepared manuscript and submit the solution by the due date on
Gradescope.

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.

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

\subsection{Problem 1: Recursive block
elimination}\label{problem-1-recursive-block-elimination}

In class, we saw ``element-wise'' variable elimination. That is, in the
first step, we eliminated the the first variable from the problem and
built a new problem for every other variable.

In fact, we can do this for multiple variables too. There's only one
small catch that has to do with pivoting, but we won't worry about that
for this problem. That is, you can assume everything you need to be
non-singular, is actually non-singular.

\textbf{Why this matters} A recent innovation are the presence of large
``tensor-cores'' that support ``matrix multiplication'' directly on GPUs
or other CPU accelerators. These can make matrix-multiplication
operations much much faster than expected. We are going to use them to
solve a linear system via block-substitution.
\href{https://github.com/corsix/amx}{As a simple example, note that
Apple's recent M-series chip includes a dedicated
``matrix-multiplication'' instruction!}
\href{https://www.intel.com/content/www/us/en/developer/articles/code-sample/advanced-matrix-extensions-intrinsics-functions.html}{As
another example, Intel has their own extensions too}

\textbf{Your problem} Write down the steps of an algorithm to solve
\(\mA \vx = \vb\) via recursive block-elimination. That is, divide
\(\mA\) into blocks: \[ \mA = \bmat{\mA_1 & \mA_2 \\ \mA_3 & \mA_4} \]
where all blocks have the same size. In this case, assume that \(\mA\)
has a power-of-two dimension so that you can do this recursively as well
until you get down to a \(1\times 1\) system. Then if
\[ \vx = \bmat{\vx_1 \\ \vx_2}, \qquad \vb = \bmat{\vb_1 \\ \vb_2} \]
are the vectors that conform to the partitioning of \(\mA\), develop an
algorithm that uses as a subroutine:

\begin{verbatim}
    # return the matrix after solving k linear systems $A^{-1} B$
    function solve(A::Matrix, B::Matrix)
    end
\end{verbatim}

and matrix-matrix multiplication to solve the linear system. Operations
you are allowed:

\begin{itemize}
\tightlist
\item
  submatrix extraction,
  e.g.~\texttt{A1\ =\ A{[}1:div(n,2),1:div(n,2){]}}
\item
  addition and subtractions of matrices or vectors,
  e.g.~\texttt{A\ .-\ B}
\item
  matrix-vector or matrix-matrix multiplication, e.g.~\texttt{A*b} or
  \texttt{A*B}
\item
  recursive calls to \texttt{solve}
\item
  you are allowed to use division only for a \texttt{n=1} case.
\end{itemize}

For testing, you can use the fact that pivoting isn't required for
symmetric positive definite matrices. So just compute a random matrix
\(\mZ\) and let \(\mA = \mZ^T \mZ\).

\textbf{Notes} You may find it helpful to think about doing this with
only \emph{two} variables at a time first. Here, I write down a few
small steps in the process to get you started. \emph{I have not fully
debugged these notes as they are meant as a supplementary guide to this
problem as a bridge from the lecture notes to the question}

Let\\
\[ \mA = \bmat{ \alpha & \beta & \vc^T \\ \gamma & \delta & \vd^T \\ \vf & \vg & \mR } \qquad \vb = \bmat{ \theta \\ \nu \\ \vh } \]
And let\\
\[ \vx = \bmat{ \omega \\ \mu \\ \vy}. \] Let
\(\mA_1 = \bmat{\alpha & \beta \\ \gamma & \delta}\). Then note that
\(\mA_1 \bmat{ \omega \\ \mu}  + \bmat{ \vc^T \\ \vd^T } \vy  = \bmat{ \theta \\ \nu }\).
So like we found in class, we can find \(\omega\) and \(\mu\) given
\(\vy\). (Assuming that \(\mA_1\) is non-singular. (Which, to be clear,
we are assuming in this problem!) Let \(\vz = \bmat{\omega \\ \mu}\).
The last expression shows that we can write \(\vz\) as a function of
\(\vy\). At this point, we can then substitute \(\vz(\vy)\) in like we
did for the other cases in the class notes.

\subsection{Problem 2: More on
Cholesky}\label{problem-2-more-on-cholesky}

\textbf{Note, we often use either \(\mF\) or \(\mL\) as the Cholesky
factor.}

Your professor also boldy asserted that we preserve positive
definiteness after one step of the Cholesky factorization. Recall that
this was: \[ \mA = \bmat{\alpha & \vb^T \\ \vb & \mC} = \mF \mF^T = 
\bmat{\gamma & 0 \\ \vf & \mF_1} \bmat{\gamma & \vf^T \\ 0 & \mF_1^T}. \]
So we set \(\gamma = \sqrt{\alpha}\) and \(\vf = \vb/\gamma\). Then, we
recursively compute the Cholesky factorization of
\[ \mC - \vb \vb^T/\alpha = \mF_1 \mF_1^T. \] This only works if
\(\mC - \vb \vb^T/\alpha\) is positive definite. This is discussed in
the notes. (Hint: you can do it from the definition of positive
definiteness: \(\vx^T \mA \vx > 0\) for \emph{all} \(\vx\).)

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  \emph{Briefly} explain why this step justifies that a Cholesky
  factorization exists for any positive definite matrix. (This is really
  like a one or two sentence answer.)
\item
  One of the most useful aspects of the Cholesky factorization is as a
  way to \emph{check} if a matrix is positive definite \emph{or}
  negative definite. A matrix is negative definite if
  \[ \vx^T \mA \vx < 0 \] for all \(\vx\). \textbf{Change} our Julia
  implementation to report if a matrix is not positive definite. (Hint:
  this relates to \emph{why} a Cholesky factorization always exists for
  a positive definite matrix in the previous problem.)
\end{enumerate}

\textbf{Use this to test} if the matrix for Poisson's equation from
Homework 1 is positive or negative definite.

\subsection{Problem 3: Direct Methods for Tridiagonal systems (Even more
on
Cholesky!)}\label{problem-3-direct-methods-for-tridiagonal-systems-even-more-on-cholesky}

In class, we said the definition of a sparse matrix was one with enough
zeros to take advantage of it. In this problem, we'll show how to take
advantage of tridiagonal structure in a matrix.

Make sure you have read through the notes on the Cholesky factorization
in our notes on Elimination methods.

Let \(\mA\) be a symmetric, tridiagonal, positive definite matrix:
\[ \mA = \bmat{\alpha_1 & \beta_1 \\
               \beta_1 & \alpha_2 & \ddots \\
                       & \ddots & \ddots & \ddots \\
                       &          & \beta_{n-2} & \alpha_{n-1} & \beta_{n-1} \\
                       &          &              & \beta_{n-1} & \alpha_n } \]

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Let \(\mA_{n-1}\) is the \(n-1 \times n-1\) leading prinicipal minor
  of \(\mA\). Suppose you are given the Choleksy factorization of
  \(\mA_{n-1} = \mL_{n-1} \mL_{n-1}^T\). Determine the structure of
  \(\mL_{n-1}\) (hint, \(\mL_{n-1} \mL_{n-1}^T\) must be tridiagonal!).
\item
  Now, show how to compute the Cholesky factorization of
  \(\mA = \mL \mL^T\) from \(\mL_{n-1}\) in a constant number of
  operations.
\item
  Use the result of the first two parts to write down an algorithm that
  computes the Cholesky factorization of \(\mA\) in a linear amount of
  work (starting from scratch!)
\end{enumerate}

\subsection{Problem 4: Direct methods for
eigenvectors.}\label{problem-4-direct-methods-for-eigenvectors.}

Suppose you try to compute an eigenvector using the ``elimination''
techniques from class. This question will ask you to investigate the
issues that arise.

\begin{verbatim}
First, we will assume we know the eigenvalue $\lambda$ and that this
\end{verbatim}

eigenvector is simple. (Don't worry if you don't know what this means.
What it means for this problem is that the linear system we build only
has the mild issue we discuss now.) Consequently, we wish to find a
solution of the linear system:
\[ \mA \vx = \lambda \vx \qquad \text{ or } \qquad (\mA - \lambda \mI) \vx = 0 . \]
This linear system is singular because any scalar multiple of an
eigenvector is also a solution. However, this might be considered a mild
singularity as we simply get to choose the scale of the solution. (If
the eigenvalue was not simple, then the solution would be more
complicated, that's why we assume simple, which formally means that
\(\rank(\mA - \lambda \mI) = n-1\))

Explain what happnes when we try and use our elimination solver
\texttt{solve1\_pivot2} using this system and also alter this routine so
that we are able to produce an eigenvector when given the matrix \(\mA\)
and \(\lambda\).

\subsection{Problem 5: QR and Cholesky}\label{problem-5-qr-and-cholesky}

Let \(\mA\) be an \(m \times n\) tall (\(m \ge n\)), full rank matrix.

Let \(\mA = \mQ \mR\) be the QR decomposition. Let \(\mF\) be the
Cholesky factor of the symmetric, positive definite matrix
\(\mB = \mA^T \mA\). (So \(\mF^T \mF = \mB\).) Show that \(\mR\) and
\(\mF\) have an extremely close relationship and discuss any uniqueness
properties of Cholesky or QR that impact this relationship.

\end{document}
