\documentclass{article}

\usepackage{fullpage}
\usepackage{amsmath,amssymb}
\usepackage{dsfont}
\usepackage{enumitem}
\usepackage{xcolor}
\usepackage{url}
\usepackage{verbatim}
\usepackage{graphicx}
\usepackage{parskip}

\usepackage[colorlinks,linkcolor=red!80!black]{hyperref}
\usepackage[noabbrev,capitalize]{cleveref}

% Use one or the other of these for displaying code.
% NOTE: If you get
%  ! Package minted Error: You must invoke LaTeX with the -shell-escape flag.
% and don't want to use minted, just comment out the next line
\usepackage{minted} \BeforeBeginEnvironment{minted}{\begingroup\color{black}} \AfterEndEnvironment{minted}{\endgroup} \setminted{autogobble,breaklines,breakanywhere,linenos}

\usepackage{listings}


% Commands for questions and answers
\definecolor{question}{rgb}{0,0,1}
\newcommand{\ask}[1]{\textcolor{question}{#1}}
\newenvironment{asking}{\begingroup\color{question}}{\endgroup}

\crefname{section}{Question}{Questions}
\newlist{qlist}{enumerate}{1}
\setlist[qlist,1]{leftmargin=*, label=\textbf{[\thesection.\arabic*]}, ref={[\thesection.\arabic*]}}
\crefname{qlisti}{Question}{Questions}

\definecolor{answer}{rgb}{0,.5,0}
\newcommand{\ans}[1]{\par\textcolor{answer}{Answer: #1}}
\newenvironment{answer}{\par\begingroup\color{answer}Answer: }{\endgroup}

\newcommand{\red}[1]{\textcolor{red}{#1}}

\definecolor{points}{rgb}{0.6,0.3,0}
\newcommand{\pts}[1]{\textcolor{points}{[#1~points]}}

\newcommand{\hint}[1]{\textcolor{black!60!white}{\emph{Hint: #1}}}
\newcommand{\meta}[1]{\textcolor{black!60!white}{\emph{#1}}}

\newcommand{\TODO}{\red{TODO}}

% misc shortcuts
\newcommand{\fig}[2]{\includegraphics[width=#1\textwidth]{#2}}
\newcommand{\centerfig}[2]{\begin{center}\includegraphics[width=#1\textwidth]{#2}\end{center}}

% Math
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\newcommand{\bigO}{\mathcal{O}}
\newcommand{\bX}{\mathbf{X}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\tp}{^\mathsf{T}}

%%begin novalidate  % overleaf code check gives false positives here
\makeatletter
% \abs{} uses auto-sized bars; \abs*{} uses normal-sized ones
\newcommand{\abs}{\@ifstar{\@Abs}{\@abs}}
\newcommand{\@abs}[1]{\left\lvert #1 \right\rvert}
\newcommand{\@Abs}[1]{\lvert #1 \rvert}

\newcommand{\norm}{\@ifstar{\@Norm}{\@norm}}
\newcommand{\@norm}[1]{\left\lVert #1 \right\rVert}
\newcommand{\@Norm}[1]{\lVert #1 \rVert}

\newcommand{\ceil}{\@ifstar{\@Ceil}{\@ceil}}
\newcommand{\@ceil}[1]{\left\lceil #1 \right\rceil}
\newcommand{\@Ceil}[1]{\lceil #1 \rceil}

\newcommand{\floor}{\@ifstar{\@Floor}{\@floor}}
\newcommand{\@floor}[1]{\left\lfloor #1 \right\rfloor}
\newcommand{\@Floor}[1]{\lfloor #1 \rfloor}

\newcommand{\inner}{\@ifstar{\@Inner}{\@inner}}
\newcommand{\@inner}[2]{\left\langle #1, #2 \right\rangle}
\newcommand{\@Inner}[2]{\langle #1, #2 \rangle}
\makeatother
%%end novalidate  % overleaf code check false-positives here


\begin{document}

\begin{center}
\Large
CPSC 440/550 Machine Learning (Jan-Apr 2024)\\
Assignment 1 --
due Friday January 19th at \textbf{5pm}
\end{center}



\textbf{IMPORTANT!!!!! Please carefully read the submission instructions that will be posted on Piazza. We will deduct up to 50\% on assignments that do not follow the instructions.}

Most of the questions below are related to topics covered in CPSC 340, or other courses listed on the prerequisite form. There are several notes available on the webpage which can help with some relevant background.

If you find this assignment to be overly difficult, that is an early warning sign that you may not be prepared to take CPSC 440
at this time. Future assignments may be longer and more difficult than this one.

We use \ask{blue} to highlight the deliverables that you must answer/do/submit with the assignment.

Cite any sources outside of course materials (including people) that you refer to,
as discussed on the syllabus slides.
Do not use ChatGPT, Copilot, or other AI assistance tools. (Something like Grammarly is okay, but not tools that suggest ``content.'')
Even if you've written code for some of these questions before (e.g.\ in 340),
\textbf{do it again}.


\clearpage
\section{Matrix Notation, Quadratics, Convexity, and Gradients \pts{20}} \label{q:mat}
\meta{Some of the notes on the course page might be useful to refresh some mathematical tools used in CPSC 340;
in particular, your instance of 340 may not have covered positive semi-definite matrices, so check those notes out if needed.}
Each part is worth \pts{2}.

\begin{qlist}
\item Consider the function
\[
f(x) = \sum_{i=1}^n\sum_{j=1}^n a_{ij}x_ix_j + \sum_{i=1}^n b_ix_i + c,
\]
where $x$ is a vector of length $n$ with elements $x_i$, $b$ is a vector of length $n$ with elements $b_i$, and $A$ is an $n \times n$ matrix with elements $a_{ij}$ (not necessarily symmetric). \ask{Write this function in matrix notation} so that it uses $A$ and $b$, and does not have summations or references to indices $i$.

\begin{answer}\TODO\end{answer}

\item \ask{Write the gradient of $f$ from the previous question, in matrix notation}.
\begin{answer}\TODO\end{answer}

\item \label{q:mat:quad-conv} \ask{Show that $f$ is convex} if $A$ is a symmetric, positive semi-definite matrix.

\begin{answer}\TODO\end{answer}

\item When $A$ is symmetric and \emph{strictly} positive definite, \ask{give a linear system whose solution minimizes $f$ in terms of $x$}.

\begin{answer}\TODO\end{answer}

\item Suppose that $A$ is symmetric and only positive \emph{semi-}definite. \ask{Will a ``good'' optimization algorithm}, e.g.\ gradient descent with an appropriate step size schedule, \ask{necessarily find one of the solutions to your linear system from the last part? If not, where else could it go?}
\hint{An example of a matrix that is psd but not strictly pd is $A = 0$.}

\begin{answer}\TODO\end{answer}

\item Consider weighted linear regression with an L2 regularizer, with regularization strength $1/\sigma^2$,
\[
f(w) = \frac12 \sum_{i=1}^n v^{(i)} (y^{(i)}  - w\tp x^{(i)} )^2 + \frac{1}{2\sigma^2}\sum_{j=1}^d w_j^2,
\]
where $\sigma > 0$ is finite and we have a vector $v$ of length $n$ containing the elements $v^i$.
As usual, $w$ is the weight vector, $x^{(i)}$ the $i$th input point, and $y^{(i)} \in \R$ its label.
\ask{Write this function in matrix notation}.

Note: you can use $\bX \in \R^{n \times d}$ as the matrix with rows $(x^{(i)})\tp$, $y$ as the vector containing the elements $y^{(i)}$, and $V$ as a diagonal matrix containing the vector $v$ along the diagonal.

\begin{answer}\TODO\end{answer}

\item Assuming that $v^i \geq 0$ for all $i$, \ask{show that $f$ from the previous part is convex}. \meta{(You can use \cref{q:mat:quad-conv} or not, as you prefer.)}

\begin{answer}\TODO\end{answer}

\item Assuming that we have $v^i \geq 0$ for all $i$, \ask{give a linear system whose solution minimizes $f$ in terms of $w$}.
\meta{(Again, use the previous results or not, your choice.)}

\begin{answer}\TODO\end{answer}


\item \label{q:mat:exp} If we fit a linear classifier with the exponential loss (used in older boosting algorithms, among other places), it would take the form
\[
f(w) = \sum_{i=1}^n \exp(-y^{(i)} w\tp x^{(i)}).
\]
\ask{Derive the gradient of this loss function.}
\begin{answer}\TODO\end{answer}


\item \label{q:mat:svr} The support vector regression objective function is
\[
f(w) = \sum_{i=1}^n \max\left\{0, \abs{w\tp x^{(i)} - y^{(i)}} - \epsilon\right\} + \frac \lambda 2 \norm{w}^2
\]
where $\epsilon$ is a non-negative hyper-parameter. It is similar to the L1 robust regression loss, but ``doesn't care'' if your prediction is less than $\epsilon$ from the target (which can reduce overfitting). \ask{Show that this non-differentiable function is convex}.

\begin{answer}\TODO\end{answer}

\end{qlist}




\clearpage
\section{Machine Learning Model Memory and Time Complexities \pts{18}}


Answer the following questions using big-O notation, and a brief explanation.
Your answers may involve $n$, $d$, and perhaps additional quantities defined in the question. 
As an example, the (linear) least squares model has $\bigO(d)$ parameters, requires $\bigO(d)$ time for prediction, and requires $\bigO(nd^2 + d^3)$ time to train.\footnote{In this course, we assume matrix operations have the ``textbook'' cost where the operations are implemented in a straightforward way, e.g.\ with \texttt{for} loops. For example, we'll assume that multiplying two $n \times n$ matrices or computing a matrix inverse simply costs $\bigO(n^3)$, rather than the $\bigO(n^\omega)$ where $\omega$ is more like $2.37$, which while true is only relevant for extremely huge matrices; $\bigO(n^3)$ is closer to the behaviour observed for actual practical matrix sizes.}

Each part is worth \pts{2}.

\begin{qlist}\color{question}
\item What is the training time for least squares (linear regression by solving the normal equations) with L2 regularization?
\begin{answer}\TODO\end{answer}

\item What is the prediction cost for a trained linear model on $t$ test examples?
\begin{answer}\TODO\end{answer}

\item What is the storage space required to make predictions with a linear regression model using Gaussian RBF features (with lengthscale, aka bandwidth, $\sigma$)?
\begin{answer}\TODO\end{answer}

\item What is the prediction time for linear regression with Gaussian RBFs (with lengthscale, aka bandwidth, $\sigma$) as features on $t$ test examples?
\begin{answer}\TODO\end{answer}

\item What is the cost of evaluating the support vector regression objective function from \cref{q:mat:svr}?
\begin{answer}\TODO\end{answer}

\item What is the cost of trying to minimize the exponential loss from \cref{q:mat:exp} by running $t$ iterations of gradient descent?
\begin{answer}\TODO\end{answer}

\item What is the cost of trying to minimize the exponential loss by running $t$ iterations of stochastic gradient descent?
\begin{answer}\TODO\end{answer}

\item What is the storage space required for the $k$-means clustering algorithm?
\begin{answer}\TODO\end{answer}

\item What is the cost of clustering $t$ examples using an already-trained k-means model?
\begin{answer}\TODO\end{answer}

\end{qlist}
\meta{Many people got very low grades on variants of this question from past years. If you're not sure how to answer questions like this, get help!}



\clearpage
\section{MAP Estimation \pts{12}}


In 340, we showed that under the assumptions of a Gaussian likelihood and Gaussian prior,
\[
y^{(i)} \sim \mathcal{N}(w^\top x^{(i)},1), \quad w_j \sim \mathcal{N}\left(0,\frac{1}{\lambda}\right),
\]
the MAP estimate is equivalent to solving the L2-regularized least squares problem
\[
f(w) = \frac{1}{2}\sum_{i=1}^n (w^\top x^{(i)} - y^{(i)})^2 + \frac \lambda 2 \sum_{j=1}^d w_j^2,
\]
in the ``loss plus regularizer'' framework.
For each of the alternate assumptions below, \ask{write it in the ``loss plus regularizer'' framework (simplifying as much as possible)}:

\begin{qlist}
\item \pts{4} Gaussian likelihood with a separate variance $\nu^{(i)} > 0$ for each training example, and Laplace prior with a separate scale $1/\lambda_j > 0$ for each variable,
\[ y^{(i)} \sim \mathcal{N}(w\tp x^{(i)},\nu^{(i)}), \quad w_j \sim \mathcal{L}\left(0,\frac{1}{\lambda_j}\right) \]
where a Laplace distribution $\mathcal L(\mu, b)$ has density $\frac{1}{2 b} \exp\left( - \abs{x - \mu} / b \right)$.

\begin{answer}\TODO\end{answer}

\item \pts{4} Robust student-$t$ likelihood and a uniform prior for $w$,
\[
p(y^{(i)} | x^{(i)}, w) =
\frac{1}{\sqrt{\nu}B\left(\frac 1 2,\frac \nu 2\right)}
\left(1 + \frac{(w\tp x^{(i)} - y^i)^2}{\nu}\right)^{-\frac{\nu + 1}{2}},
\quad w \sim \operatorname{Uniform}\left( \{ w : \norm w \le C \} \right),
\]
where $B$ is the ``Beta" function, and the parameter $\nu$ is called the ``degrees of freedom.''\footnote{This likelihood is more robust than the Laplace likelihood, but leads to a non-convex objective.}
You can use $V_d$ as the volume of the $d$-dimensional unit ball $\{ w : \norm w \le 1 \}$.
\hint{Think carefully about what a uniform prior density means for the posterior density of $w$.}

\begin{answer}\TODO\end{answer}

\item \pts{4} If $y^{(i)}$ represents counts, we can use a Poisson-distributed likelihood, and take a Gaussian prior on $w$ centered at $u \in \R^d$:
\[
p(y^{(i)} | x^{(i)}, w) = \frac{\exp(y^{(i)} w\tp x^{(i)})\exp(-\exp(w^\top x^{(i)}))}{y^{(i)}!}, \quad
w_j \sim \mathcal N\left(u_j, \frac1\lambda\right)
\]
\begin{answer}\TODO\end{answer}

\end{qlist}



\clearpage
\section{K-Means Clustering \pts{15}}

\meta{This and following questions use code available from the course webpage; get it from \texttt{a1.zip}.}

\meta{You'll need Python 3 and the packages \texttt{numpy}, \texttt{scipy}, and \texttt{matplotlib}.
If you don't have them installed already,
install them with \texttt{conda install}, your system package manager, or
\texttt{pip install numpy scipy matplotlib}.}
 

If you run \texttt{main.py kmeans}, it will load a dataset with two features and a very obvious clustering structure. It will then apply the $k$-means algorithm with a random initialization. The result of applying the algorithm will thus depend on the randomization, but a typical run might look like this:\\
\centerfig{.5}{figs/kmeans-bad.png}
(Note that the colours are arbitrary due to the label switching problem.)
But the ``correct'' clustering (that was used to make the data) is something more like this:\\
\centerfig{.5}{figs/kmeans-good.png}

If you run the demo several times, it will find different clusterings. To select among clusterings for a \emph{fixed} value of $k$, one strategy is to minimize the sum of squared distances between examples $x^{(i)}$ and their means $w_{y^{(i)}}$,
\[
f(w^{(1)},w^{(2)},\dots,w^{(k)},y^{(1)},y^{(2)},\dots,y^{(n)}) = \sum_{i=1}^n \norm{x^{(i)} - w^{( y^{(i)} )}}_2^2 = \sum_{i=1}^n \sum_{j=1}^d \left(x^{(i)}_j -  w^{( y^{(i)} )}_j\right)^2.
\]
where $y^{(i)} \in \argmin_{c \in [n]} \norm{x^{(i)} - w^{(c)}}$ is the index of the closest mean to $x_i$. This is a natural criterion because the steps of k-means alternately optimize this objective function in terms of the $w^{(c)}$ and the $y^{(i)}$ values.
 
\begin{qlist}
 \item \pts{3} Complete the function \texttt{KMeans.loss}, inside \texttt{k\_means.py}, that takes in a data matrix $\bX$, a vector of corresponding cluster assignments $y$, and a matrix of cluster means $W$, and computes this objective function. \ask{Hand in your code.}
\begin{answer}\TODO\end{answer}


 \item \pts{1} Modify your code to, instead of/in addition to printing the number of labels that change on each iteration, print the value of \texttt{KMeans.loss} after each iteration. \ask{What trend do you observe?} (No need to hand in code or specific loss values for this, just describe the trend.)
\begin{answer}\TODO\end{answer}

 \item \pts{2} \texttt{main.py kmeans-best} will call the function \verb|q_kmeans_best()| in \texttt{main.py}, which calls the \verb|best_fit| function to rerun $k$-means 50 times and take the one with the lowest error. Complete the \verb|best_fit| function, and \ask{hand in your code and the resulting plot}.
\begin{answer}\TODO\end{answer}

 \item \pts{3} \ask{Explain why} the \texttt{KMeans.loss} function should not be used to choose $k$ -- even if you evaluate it on test data.
\begin{answer}\TODO\end{answer}


 \item \pts{1} The data in \texttt{cluster\_data\_2.npz} is exactly the same as \texttt{cluster\_data.npz}, except that it has four outliers that are far away from the data. \texttt{main.py kmeans-outliers} will run your code from above to find the best of 50 runs on this data and save it as \texttt{figs/kmeans-outliers.png}. \ask{Hand in this plot; are you satisfied with it?}
 \begin{answer}\TODO\end{answer}

 \item \pts{5} Implement the $k$-\emph{medians} algorithm in \texttt{kmedians.py}, which assigns examples to the nearest $w^{(c)}$ in the L1-norm, and then updates all the $w^{(c)}$ by setting them to the ``median'' of the points assigned to the cluster (defining the $d$-dimensional median as the concatenation of the medians along each dimension). For this algorithm it makes sense to use the L1-norm version of the error (where $y^{(i)}$ now represents the closest centre in the L1-norm),
\[
f(w^{(1)},w^{(2)},\dots,w^{(k)},y^{(1)},y^{(2)},\dots,y^{(n)}) = \sum_{i=1}^n \norm{x^{(i)} - w^{(y^{(i)})}}_1 = \sum_{i=1}^n \sum_{j=1}^d |x^{(i)}_j - w^{(y^{(i)})}_j|,
\] 
 \texttt{main.py kmedians-outliers} will find the best of 50 models with $k = 4$ for you, once you've finished the \texttt{KMedians} class.
 \ask{Hand in your code and plot}. \ask{Is this better?}
 \begin{answer}\TODO\end{answer}
\end{qlist}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Regularization and Hyper-Parameter Tuning \pts{20}}

If you run \verb|main.py lsq|, it will:
\begin{enumerate}
\item Load a one-dimensional regression dataset.
\item Fit a least-squares linear regression model.
\item Draw a figure showing the training/testing data and what the model looks like.
\end{enumerate}
Unfortunately, this is not a great model of the data, and the figure shows that a linear model with a $y$ intercept of 0 is probably not suitable.

\begin{qlist}
\item \label{q:lsq:plain} \pts{5} Finish the \verb|LeastSquaresBias| class, in \texttt{least\_squares.py}.
Using this class should fit and predict with the model
\[
y^{(i)} = w\tp x^{(i)} + b
,\]
where both $w$ and $b$ are fit to data with least squares;
implement it however you like (there are a few reasonable options).
\texttt{main.py lsq-bias} will then run it for you.
\ask{Hand in your new class and the updated plot.}
\begin{answer}\TODO\end{answer}

\item \pts{6} Allowing a non-zero y-intercept improves the prediction substantially, but the model still seems sub-optimal because there are obvious non-linear effects. Complete the model \verb|leastSquaresRBFL2| that implements \emph{least squares using Gaussian radial basis functions (RBFs) with L2-regularization}.
\texttt{main.py lsq-rbf} will then run it for you.

Use \texttt{lam} for the L2 regularization parameter\footnote{We can't name it \texttt{lambda}, because that's a reserved word in Python. You could technically type $\lambda$, but there are a lot of Unicode characters for $\lambda$ that all look really similar, so it gets annoying\dots.}, and \texttt{sigma} for the lengthscale of the Gaussian features. Note that your L2 regularization should correspond to minimizing the loss function $\norm{\bX w - y}^2 + \lambda \norm{w}^2$; some versions instead correspond to putting a $\frac1n$ in front of the loss term.

\ask{Hand in your function and the plot generated with $\lambda = 1$ and $\sigma = 1$.}

\hint{The function \texttt{utils.euclidean\_dist\_squared}, which was also used in our k-means implementation, efficiently computes the squared Euclidean distance between all pairs of rows in two matrices.}

\begin{answer}\TODO\end{answer}

\item \pts{6} \label{q:lsq:split}
Modify the script, in the function \verb|lsq-rbf-split|, to split the training data into a ``train'' and ``validation'' set (you can use half the examples for training and half for validation), and use these to select $\lambda$ and $\sigma$ from some reasonable set of candidate values. You'll probably want to vary these by a few orders of magnitude either smaller or larger from $1$.

(Although in real work you'd probably use some pre-coded helpers for this, code it yourself here.)

\ask{Hand in your modified function, the selected $\lambda$ and $\sigma$, and the plot you obtain by refitting on the full dataset with the best values of $\lambda$ and $\sigma$.}

\begin{answer}\TODO\end{answer}


\item \pts{3} Consider (but you don't have to implement) a model combining the first two parts of this question,
\[
y^{(i)} = w\tp x^{(i)} + b + v\tp z^{(i)},
\]
where $z^{(i)}$ are the Gaussian RBFs and $v$ is a vector of regression weights for those features. Suppose that we first fit $w$ and $b$ assuming that $v=0$ as in \cref{q:lsq:plain}, and then we fit $v$ with $w$ and $w_0$ fixed (you could use your code for \cref{q:lsq:split} to do this by modifying the $y^{(i)}$ values). \ask{Why would someone want to do this?}

\hint{Think about how this model would behave on a test point $\tilde x = 10$.}

\begin{answer}\TODO\end{answer}
\end{qlist}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Optimal matrix multiplication \pts{15}}

Suppose we have a bunch of matrices $M_1, M_2, \dots, M_k$,
where $M_i$ is of shape $d_{i-1} \times d_{i}$,
and we want to find the $d_0 \times d_{k}$ matrix $M_1 M_2 \cdots M_k$.
There are lots of ways to compute this, since matrix multiplication is associative:
for example, if $k=6$
we could do $M_1 (M_2 (M_3 (M_4 (M_5 M_6))))$,
$M_1 ((M_2 M_3) ((M_4 M_5) M_6))$,
etc. The number of possibilities is exponential in $k$.

Recall that
if $A$ is $m \times n$ and $B$ is $n \times p$,
we can compute $A B$ with $\bigO(m n p)$ scalar operations;
this is typically a ``good-enough'' model for how much this costs in practice for the kinds of problems that appear in ML.

    Fill in the \texttt{mmul\_order} function in \texttt{mmul\_order.py}
    to use dynamic programming to find the optimal order of multiplications.

    Make sure you follow the comments there as to the input and output format of your function,
    and that your algorithm is $\bigO(k^3)$ or better.

    \texttt{main.py mmul-order} will call the function for you with a few test cases;
    \ask{hand in your code and the cost and order for the final part}.
\begin{answer}\TODO\end{answer}








\end{document}
