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

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



\begin{document}
\maketitle

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

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.

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

** 2023-09-20: Question 5 was updated from Gradient Descent to Steepest
Descent **

\hypertarget{week-4}{%
\section{Week 4}\label{week-4}}

Plan to do these problems during the first week.

\hypertarget{problem-1-getting-poissons-equations-2d-mesh-to-converge-with-richardson}{%
\subsection{Problem 1: Getting Poisson's equations (2d mesh) to converge
with
Richardson}\label{problem-1-getting-poissons-equations-2d-mesh-to-converge-with-richardson}}

For this problem, we will reuse the matrix for Poisson's equation from
HW1. Please use the version of the 2d mesh with the boundary conditions
eliminated. That is, you should only have unknowns for the nodes on the
interior. Also, make sure you solve a system where the diagonal elements
are all positive so that it is symmetric positive definite.

(This may require flipping the sign of the system.)

If you have questions about how to get this right, please do ask during
office hours with specific questions.

In class we showed that the Richardson method for \(\mA \vx = \vb\)
\[ \vx\itn{k+1} = \vx\itn{k} + \alpha \vr\itn{k} \] will converge for
any positive semi definite matrix as long as \(\alpha < 2 / \rho(\mA)\).
We have not yet determined ways of computing \(\rho(\mA)\), however.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Provide evidence that setting \(\alpha = 1\) does not result in a
  method that converges for Poisson's equation with \(n = 10\).
\item
  Experimentally or theoretically, find a value of \(\alpha\) such that
  using the Richardson method will converge for Poisson's equation with
  \(n=10\).
\item
  What happens when you try \(n=20\), does the same value of \(\alpha\)
  work? If not, find a new one that does.
\item
  Experimentally, find a value of \(\alpha\) that takes the fewest
  iterations to produce a solution with relative residual \(10^{-5}\)
  for both \(n=10\) and \(n=20\).
\item
  Prove that using \(\alpha < 1/4\) will result in a method that
  converges.
\end{enumerate}

\hypertarget{problem-2-matrix-and-vector-norms}{%
\subsection{Problem 2: Matrix and vector
norms}\label{problem-2-matrix-and-vector-norms}}

Consider the following function: \[ f(\mA) = \max_{i,j} |A_{i,j}|. \]

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Show that \(f\) is a matrix norm. (Very easy!)
\item
  Show that \(f\) does not satisfy the sub-multiplicative property.
\item
  Show that there exists \(\sigma > 0\) such that:
  \[ g(\mA) = \sigma f(\mA) \] is a sub-multiplicative matrix-norm.
\end{enumerate}

\hypertarget{problem-3-solving-the-pagerank-linear-system-of-equations-via-richardson.}{%
\subsection{Problem 3: Solving the PageRank linear system of equations
via
Richardson.}\label{problem-3-solving-the-pagerank-linear-system-of-equations-via-richardson.}}

We will derive the following at some point in the class, but you have
everything you need to solve this right now. Let \(\mP\) be a matrix
that is entry-wise non-negative (\(P_{i,j} \ge 0\)) and also has columns
that sum to \(1\) (\(\sum_{i} P_{i,j} = 1\)).

Informally, the matrix \(\mP\) describes the probability of moving
around a graph, just like the matrix \(\mT\) described the probability
of moving around the game of Chutes and Ladders. The interpretation of
\(P_{i,j}\) is the probability of moving from node \(j\) to node \(i\).

The PageRank linear system is:
\[ (\mI - \alpha \mP) \vx = (1-\alpha) \vv \] where \(v_i \ge 0\) and
\(\sum_i v_i = 1\). The solution \(x_{i}\) can be interpreted as the
asymptotic fraction of time a random walk spends at a node of a graph
when, at each step, it moves in the graph (according to \(\mP\)) with
probability \(\alpha\) and it \emph{jumps} according to \(\vv\) with
probability \((1-\alpha)\).

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Show that we can solve the PageRank linear system with
  \(0 < \alpha < 1\) using the Richardson method.
\item
  Use your sparse matrix tools from HW1 or other tools you implement
  yourself to explain how to implement the Richardson method for
  PageRank on the graph provided in files. Your code should allow us to
  specify:

  \begin{itemize}
  \tightlist
  \item
    value of alpha (for PageRank) and value of \(\alpha\) for Richardson
  \item
    relative residual
  \item
    maximum number of iterations
  \end{itemize}

  You should be able to run this code on large sparse matrices with
  hundreds of thousands of nodes and millions of edges. (As you will do
  in the next step.)

  Your code must execute directly on CSC or CSR arrays and cannot use
  any built in Julia functions for sparse matrices.
\item
  How many iterations does it take to compute a relative residual of
  less than \(10^{-5}\) on the data in files with \(\alpha = 0.85\) (for
  PageRank) and \(\alpha = 0.5\) and \(\alpha = 1.0\) for Richardson.

  Simple random test data

\begin{verbatim}
 using SparseArrays, Random
 Random.seed!(0)
 P = sprand(1_000_000,1_000_000, 15/1_000_000) |> 
         A->begin fill!(A.nzval,1); A; end |> 
         A->begin ei,ej,ev = findnz(A); d = sum(A;dims=2); 
                 return sparse(ej,ei,ev./d[ei], size(A)...); end
 colptr,rowval,nzval = P.colptr, P.rowval, P.nzval
\end{verbatim}

  Actual data to report

\begin{verbatim}
 using ZipFile, DelimitedFiles,SparseArrays
 function load_data()
   r = ZipFile.Reader("wikipedia-2005.zip")
   try
     @assert length(r.files) == 1
     f = first(r.files)
     data = readdlm(f,'\n',Int)
     n = data[1]
     colptr = data[2:2+n] # colptr has length n+1
     rowval = data[3+n:end] # n+2 elements before start of rowval
     A = SparseMatrixCSC(n,n,colptr,rowval,ones(length(rowval))) |> 
       A->begin ei,ej,ev = findnz(A); d = sum(A;dims=2); 
       return sparse(ej,ei,ev./d[ei], size(A)...) end
   finally
     close(r)
   end
 end
 P = load_data()
 colptr,rowval,nzval = P.colptr, P.rowval, P.nzval
\end{verbatim}

  Use the uniform vector \(\vv = \ve/n\) where \(n\) is the row
  dimension of the matrix \(\mP\).

  On the file
  \url{https://www.cs.purdue.edu/homes/dgleich/cs515-2023/homeworks/wikipedia-2005.zip}
\end{enumerate}

\hypertarget{problem-4-back-to-chutes-and-ladders.}{%
\subsection{Problem 4: Back to Chutes and
Ladders.}\label{problem-4-back-to-chutes-and-ladders.}}

For this problem, we'll use the linear system for Chutes and Ladders
from HW1, 6-1.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Develop a report on using Richardson to solve the linear system for
  the Chutes and Ladders game. Relevant dimensions include things that
  we have discussed in previous problems.

  \begin{itemize}
  \tightlist
  \item
    Convergence
  \item
    Overall speed
  \item
    Choices of alpha
  \end{itemize}
\item
  \textbf{More challenging} In our HW1 Chutes and Ladders system, we had
  two ways of computing the expected length of the game (HW1 6-1 and HW1
  6-2). These should have given you the same result. (They do!) These
  were different as one involved \(\mT\) and the other involved
  \(\mT^T\). Use the Neumann series to show that these are actually the
  same (mathematically).
\end{enumerate}

\hypertarget{week-5}{%
\section{Week 5}\label{week-5}}

Plan to do these problems the second week.

\hypertarget{problem-5-implementing-steepest-descent}{%
\subsection{Problem 5: Implementing steepest
descent}\label{problem-5-implementing-steepest-descent}}

For steepest descent, the most expensive computational step is computing
matrix-vector products with the matrix \(\mA\). Show how to implement
this algorithm with only a single matrix-vector product per iteration.
Your algorithm should stop when the residual is sufficiently small. Your
algorithm, including checking the stopping condition, can only use one
matrix-vector product per iteration.

(Hint: you will have to figure out how to \emph{update} the expression
for the gradient and solution and reformulate the algorithm!)

\hypertarget{problem-6-implementing-coordinate-descent}{%
\subsection{Problem 6: Implementing coordinate
descent}\label{problem-6-implementing-coordinate-descent}}

For coordinate descent, the most expensive step is computing
\([\vg_k]_i\). Use the routines -- or small modifications -- from
Homework 1 to compute this term efficiently by using the compressed
sparse row or compressed sparse column arrays. Develop an algorithm to
solve \(\mA \vx = \vb\) via coordinate descent with two strategies to
pick the coordinates: (1) by cycling through 1 to \(n\) (called cyclic
coordinate descent) and (2) by randomly picking a coordinate (called
random coordinate descent). If row (or column \(i\)) of the matrix has
\(k\) non-zeros, your algorithm must do work proportional to \(k\) in
each iteration.

Your solution should justify your implementation is correct.

\hypertarget{problem-7-generalizing-coordinate-and-gradient-descent}{%
\subsection{Problem 7: Generalizing coordinate and gradient
descent}\label{problem-7-generalizing-coordinate-and-gradient-descent}}

The algorithms for coordinate and gradient descent are for solving
\(\mA \vx = \vb\) on a symmetric positive definite matrix \(\mA\)
correspond to trying to minimize
\[ f(\vx) = \tfrac{1}{2} \vx^T \mA \vx - \vx^T \vb \] by, at each step,
exactly minizing this function for a single coordinate (coordinate
descent) or by exactly this function in the direction of the gradient
(gradient descent). Let \(\vx_k\) represent the approximate solution
after \(k\) steps. Let \(\vg_k\) represent the gradient
\(\mA \vx_k - \vb\) at the \(k\)th step. Then the algorithms can be
written:
\[ \vx_{k+1} = \vx_{k} - \frac{[\vg_k]_i}{A_{i,i}} \ve_i \qquad \vx_{k+1} = \vx_{k} - \frac{\vg_k^T \vg_k}{\vg_k^T \mA \vg_k} \vg_k \]
respectively.

Suppose that we wish to generalize coordinate descent to \emph{two}
coordinates. Then the update can be written:
\[ \vx_{k+1} = \vx_k + \mE_{ij} \vc_k \] where
\(\mE_{ij} = \bmat{ \ve_i & \ve_j }\) and \(\vc_k = \bmat{c_i \\ c_j}\).
That is, we have to choose two coefficients: \(c_i, c_j\) in order to
figure out the next step. Like coordinate descent, we are going to pick
these to exactly minimize \(f(\vx_{k+1})\) as a function of
\(c_i, c_j\).

Give the resulting iterative algorithm.

\hypertarget{problem-8-reporting-on-algorithms}{%
\subsection{Problem 8: Reporting on
algorithms}\label{problem-8-reporting-on-algorithms}}

Use your software from problem 5 and 6 to prepare reports analyzing the
behavior, convergence, and work on the following linear systems:

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  The 2d Laplacian system from Problem 1 on this homework with \(n=40\)
\item
  The Chutes and Ladders linear system. (Note that this is not symmetric
  positive definite, so the report may be different.)
\end{enumerate}

compare the performance of the algorithms from problems 5 and 6 in terms
of the \emph{norm of the relative residual} compared with the
\emph{number of non-zeros used}. Use a relative residual of \(10^{-4}\).

We are not comparing by number of iterations because each iteration does
vastly different amounts of work. An iteration of coordinate descent
does a very small amount of work, and we need \(n\) iterations of cyclic
coordinate descent to give us the same amount of work as 1 iteration of
gradient descent. Hence, we'll use the total number of non-zero entries
used. So if a row/column used in coordinate descent has \(k\) non-zero
entries, then we increase our work count by \(k\). For each iteration of
gradient descent, we use the total number of non-zeros in \(\mA\) work.

Your answer should probably have both a table and a plot of how the
residual descreases with \emph{work}. You should show how the residual
decreases from \(10^{-1}\) to \(10^{-4}\) as a function of work. This
suggests a log-scaling of a \(y\)-axis.

\hypertarget{problem-9-viral-spreading-and-eigenvalues.}{%
\subsection{Problem 9: Viral spreading and
eigenvalues.}\label{problem-9-viral-spreading-and-eigenvalues.}}

Use what you know about eigenvalues to help address or answer the
prompt:

\begin{quote}
For viral spreading with the linear algebraic approximation
(\texttt{approx\_evolve\_steps}), the largest eigenvalue of \(\rho \mA\)
is relevant to the problem because
\end{quote}

No more than a few sentences and perhaps a small figure is needed.

\hypertarget{problem-10-jacobi-norms-and-diagonal-dominance}{%
\subsection{Problem 10: Jacobi, Norms, and Diagonal
Dominance}\label{problem-10-jacobi-norms-and-diagonal-dominance}}

One of the common classes of matrices are called diagonally dominant.
There are a number of subclasses of diagonal dominance. Here, we will
study a related property. Consider a general linear system of equations
\[ \mA \vx = \vb \] Let \(\mA = \mD + \mN\) be a diagonal and
off-diagonal splitting of \(\mA\). In this case
\[ \mD = \text{ diagonal entries of } \mA \]
\[ \mN = \mA - \mD = \text{ all non-diagonal entries of } \mA \] Suppose
that we know that \[ \normof{\mN \vx} \le \gamma \normof{\mD \vx} \] for
the vector 2-norm and \(\gamma < 1\) for all vectors \(\vx\). Show that
the Jacobi method converges. (This is diagonally dominant in the sense
that the the diagonal terms dominant the off-diagonal terms in a
2-norm.)

\end{document}
