\documentclass{article} \usepackage[utf8]{inputenc} \usepackage{amsmath,amssymb} \usepackage{algorithm,algpseudocode} \usepackage{bm} \usepackage{hyperref} \usepackage{float} \usepackage{subcaption} \usepackage{enumitem} \usepackage{xcolor} \definecolor{darkgreen}{rgb}{0.0, 0.5, 0.0} \definecolor{aaltoBlue}{RGB}{0, 51, 102} \definecolor{pinegreen}{cmyk}{0.92,0,0.59,0.25} \definecolor{royalblue}{cmyk}{1,0.50,0,0} \definecolor{lavander}{cmyk}{0,0.48,0,0} \definecolor{violet}{cmyk}{0.79,0.88,0,0} \definecolor{lightblue}{RGB}{173, 216, 230} \usepackage{tikz} \usepackage{pgfplots} \usetikzlibrary{calc,arrows.meta,positioning,fit,shapes.geometric,shapes.misc,intersections,patterns,decorations.pathreplacing,backgrounds,matrix} \usepgfplotslibrary{fillbetween,groupplots} \setcounter{secnumdepth}{-\maxdimen} % Long unbreakable tokens (the ISO/IEC author names in citations, term names % in the title) otherwise overflow the margin by a few points; a modest % emergency stretch lets TeX loosen the paragraph instead. \emergencystretch=2em \usepackage{mleftright} % Project macros used by this entry (auto-extracted). \newcommand{\bayeshypothesis}{h^{*}} \newcommand{\condnumber}[1]{\kappa\!\left(#1\right)} \newcommand{\covmtx}[1]{\mathbf{C}^{(#1)}} \newcommand{\covvec}[1]{\mathbf{c}^{(#1)}} \newcommand{\dataset}{\mathcal{D}} \newcommand{\eigval}[1]{\lambda_{#1}} \newcommand{\expect}{\mathbb{E} } \newcommand{\feature}{x} \newcommand{\featuredim}{d} \newcommand{\featureidx}{j} \newcommand{\featuremtx}{\mX} \newcommand{\featurevec}{\vx} \newcommand{\fixedpointop}{\mathcal{F}} \newcommand{\gdstep}[1]{\mathcal{F}^{(#1)}} \newcommand{\hypothesis}{h} \newcommand{\iteridx}{t} \newcommand{\labelvec}{\vy} \newcommand{\lrate}{\eta} \newcommand{\mI}{\mathbf{I}} \newcommand{\mM}[0]{{\bf M}} \newcommand{\mX}[0]{{\bf X}} \newcommand{\normgeneric}[2]{\mleft\lVert #1 \mright\rVert_{#2}} \newcommand{\nrfeatures}{\featuredim} \newcommand{\pair}[2]{\left( #1,#2 \right)} \newcommand{\perturbation}[1]{{\bm \varepsilon}^{(#1)}} \newcommand{\pinv}[1]{{#1}^{+}} \newcommand{\probdist}{\mathbb{P}} \newcommand{\reals}{\mathbb{R}} \newcommand{\regparam}{\alpha} \newcommand{\sampleidx}{r} \newcommand{\samplesize}{m} \newcommand{\trainset}{\dataset^{(\rm train)}} \newcommand{\truelabel}{y} \newcommand{\vb}[0]{{\bf b}} \newcommand{\vw}[0]{{\bf w}} \newcommand{\vx}[0]{{\bf x}} \newcommand{\vy}[0]{{\bf y}} \newcommand{\weights}{\vw} \begin{document} \section{linear regression} \subsection*{Author and editors} Alexander Jung\textsuperscript{*}, Department of Computer Science, Aalto University, Espoo, Finland (ORCID \mbox{\href{https://orcid.org/0000-0001-7538-0990}{0000-0001-7538-0990}}) \\ Editors: Konstantina Olioumtsevits (ORCID \mbox{\href{https://orcid.org/0000-0003-2697-8659}{0000-0003-2697-8659}}), Ekkehard Schnoor (ORCID \mbox{\href{https://orcid.org/0000-0001-8252-1499}{0000-0001-8252-1499}}), Aalto University, Espoo, Finland \\ \textsuperscript{*}Corresponding author: \href{mailto:alex.jung@aalto.fi}{alex.jung@aalto.fi} \subsection{Abstract} Linear regression is a regression method that learns a linear hypothesis map for predicting the numeric label of a data point from its features. Its least squares variant chooses a map that minimizes the average squared error loss on a training set. Linear regression is an instance of empirical risk minimization (ERM) that is obtained by using the linear model and the squared error loss. The optimal model parameters are determined by the normal equations. This system can be solved in closed form or iteratively, for example, by gradient descent (GD). Adding a penalty term to the average squared error loss yields regularized variants such as ridge regression and least absolute shrinkage and selection operator (Lasso). \subsection{Definition} Linear regression can be used to predict tomorrow's temperature from weather measurements recorded over the last five days, using them as features and tomorrow's temperature as the label. Linear regression methods use a linear model: they learn a linear hypothesis map that delivers a prediction of the numeric label of a data point. The prediction is based solely on the features of the data point, which are fed as input to the learned hypothesis map. \par Formally, linear regression learns a linear hypothesis map $\hypothesis^{(\weights)}(\featurevec) = \weights^{\top}\featurevec$ to predict the numeric label $\truelabel \in \reals$ of a data point from its feature vector $\featurevec = \big(\feature_{1},\,\ldots,\,\feature_{\nrfeatures}\big)^{\top} \in \reals^{\nrfeatures}$. Here, $\nrfeatures$ denotes the number of features of a data point. The model parameters $\weights$ are learned from a training set $\trainset = \big\{ \pair{\featurevec^{(\sampleidx)}}{\truelabel^{(\sampleidx)}} \big\}_{\sampleidx=1}^{\samplesize}$ of data points. Appending a constant feature to $\featurevec$ lets one entry of $\weights$ act as an intercept, so $\weights^{\top}\featurevec$ represents an affine function of the original measurements. Below, $\featurevec$ denotes whichever feature vector is fed to the linear hypothesis map, original or augmented, unless stated otherwise. \par The least squares variant of linear regression measures the quality of a linear hypothesis map by its average loss, the squared error loss, on the training set. As an instance of empirical risk minimization (ERM), it learns the model parameters $\weights$ by solving the optimization problem \[ \min_{\weights \in \reals^{\nrfeatures}} \frac{1}{\samplesize} \sum_{\sampleidx=1}^{\samplesize} \big( \truelabel^{(\sampleidx)} - \weights^{\top} \featurevec^{(\sampleidx)} \big)^{2} \text{.} \] Fig.~\ref{fig_linreg_dict} illustrates this linear regression problem for a training set with a constant scalar feature $\feature^{(\sampleidx)}=1$ for $\sampleidx=1,\ldots,\samplesize$. \begin{figure}[H] \centering \begin{tikzpicture}[scale=0.7, y=1.1cm, x=0.5cm] \foreach \x/\y in { 1/2, 4/3, 7/4 } { \draw[dashed, gray] (\x, 3) -- (\x, \y); \filldraw[blue] (\x, \y) circle (2pt); } \draw[dashed, thick] (0.5, 3) -- (10.5, 3) node[right] {$\widehat{w}=3$}; \end{tikzpicture} \caption{For a linear model with $\nrfeatures=1$ and the constant feature $\feature=1$ for every data point (blue filled circles), linear regression reduces to computing the average $\widehat{w} = (1/\samplesize) \sum_{\sampleidx=1}^{\samplesize} \truelabel^{(\sampleidx)}$ of the labels: this average minimizes the average squared error loss. The horizontal dashed line sits at the height $\widehat{w}$, and the vertical gray segments are the errors $\truelabel^{(\sampleidx)} - \widehat{w}$ entering that loss \label{fig_linreg_dict}} \end{figure} The optimization problem can be written more compactly using the feature matrix $\featuremtx$ and the label vector $\labelvec$, \[ \featuremtx = \big(\featurevec^{(1)},\,\ldots,\,\featurevec^{(\samplesize)}\big)^{\top} \in \reals^{\samplesize \times \nrfeatures} \text{,} \qquad \labelvec = \big( \truelabel^{(1)},\,\ldots,\,\truelabel^{(\samplesize)} \big)^{\top} \in \reals^{\samplesize} \text{.} \] In terms of $\featuremtx$ and $\labelvec$, the optimization problem reads \[ \min_{\weights \in \reals^{\nrfeatures}} \underbrace{\frac{1}{\samplesize} \normgeneric{\labelvec - \featuremtx \weights}{2}^{2}}_{f(\weights)} \text{.} \] The average squared error loss $f(\weights)$ is a convex and differentiable function of $\weights$. By the zero-gradient condition (Boyd and Vandenberghe, 2004, Sect.~4.2.3), a vector $\widehat{\weights}$ solves this optimization problem if and only if it satisfies the system of linear equations (Strang, 2016) \begin{equation} \label{eq_linreg_normal_eq_dict} \featuremtx^{\top}\featuremtx \widehat{\weights} = \featuremtx^{\top} \labelvec \text{.} \end{equation} The normal equations \eqref{eq_linreg_normal_eq_dict} always have a solution because $\featuremtx^{\top}\labelvec$ lies in the column space of $\featuremtx^{\top}$, which coincides with the column space of $\featuremtx^{\top}\featuremtx$ (Meyer, 2000, Sect.~4.5). \par However, the solution of the normal equations \eqref{eq_linreg_normal_eq_dict} is unique only if the feature matrix $\featuremtx$ has full column rank. This requires at least as many data points as features, $\samplesize \geq \nrfeatures$. When this holds, and the columns of $\featuremtx$ are linearly independent, the matrix $\featuremtx^{\top}\featuremtx$ is invertible and the normal equations \eqref{eq_linreg_normal_eq_dict} have the unique closed-form solution $\widehat{\weights} = \big(\featuremtx^{\top}\featuremtx\big)^{-1} \featuremtx^{\top} \labelvec$ (Golub and Van Loan, 2013, Sect.~5.3.1). \par If instead $\samplesize < \nrfeatures$ (more features than data points), $\featuremtx$ cannot have full column rank, $\featuremtx^{\top}\featuremtx$ is singular, and the optimization problem has infinitely many solutions (Golub and Van Loan, 2013, Sect.~5.5). All of them incur the same minimal average squared error loss on the training set, but their predictions for data points outside the training set can differ arbitrarily. Choosing among them is therefore a matter of generalization (see overfitting). A unique solution can be selected, for example, by picking the minimum-norm solution $\pinv{\featuremtx} \labelvec$ given by the pseudoinverse $\pinv{\featuremtx}$ (Golub and Van Loan, 2013, Sect.~5.5.2). \par A common remedy is regularization: adding a penalty term to the average squared error loss. The penalty term can be interpreted as an estimate of how much higher the loss is on data points outside the training set than on it (Hastie et al., 2009, Sect.~7.4). This yields regularized variants of linear regression. Ridge regression uses the penalty term $\regparam \normgeneric{\weights}{2}^{2}$ and the least absolute shrinkage and selection operator (Lasso) uses $\regparam \normgeneric{\weights}{1}$. The hyperparameter $\regparam > 0$ sets the weight of the penalty term. \par Linear regression also has a statistical interpretation. Consider a probabilistic model with a joint probability distribution $\probdist^{(\featurevec, \truelabel)}$ over the features and the label. Under the squared error loss, the Bayes estimator (i.e., the hypothesis with minimum risk) is the conditional expectation $\expect\{\truelabel \mid \featurevec\}$ (Lehmann and Casella, 1998, Ch.~4; Papoulis and Pillai, 2002). In this statistical interpretation, $\featurevec$ denotes the original feature vector. Let the features and the label be jointly Gaussian random variables (Gaussian RVs) with zero mean and invertible covariance matrix $\covmtx{\featurevec} = \expect\{\featurevec \featurevec^{\top}\}$. This Bayes estimator is then (Bishop, 2006, Sect.~2.3.1, Eq.~(2.81)) \[ \bayeshypothesis(\featurevec) = \big(\weights^{\star}\big)^{\top}\featurevec \text{, with } \weights^{\star} = \big(\covmtx{\featurevec}\big)^{-1}\, \covvec{\featurevec,\truelabel} \text{.} \] Here, $\covvec{\featurevec,\truelabel} = \expect\{\featurevec \truelabel\}$ is the covariance between the features $\featurevec$ and the label $\truelabel$. This Bayes estimator is linear in $\featurevec$. Without the zero-mean assumption, it is an affine function of $\featurevec$, which a constant feature absorbs as before. \par Least-squares linear regression is the ERM counterpart of this predictor. It replaces $\covmtx{\featurevec}$ and $\covvec{\featurevec,\truelabel}$ with the sample-based estimates $(1/\samplesize)\featuremtx^{\top}\featuremtx$ and $(1/\samplesize)\featuremtx^{\top}\labelvec$. For invertible $\featuremtx^{\top}\featuremtx$, this recovers the unique solution $\widehat{\weights} = \big(\featuremtx^{\top}\featuremtx\big)^{-1}\featuremtx^{\top}\labelvec$ of the normal equations \eqref{eq_linreg_normal_eq_dict}. \par Instead of solving the normal equations \eqref{eq_linreg_normal_eq_dict} directly (via the inverse or pseudoinverse), machine learning (ML) methods often solve it by an iterative optimization method, as implemented in widely used ML software libraries such as scikit-learn (Pedregosa et al., 2011) and PyTorch (Paszke et al., 2019). Starting from initial model parameters $\weights^{(0)}$, such a method repeatedly applies an update operator $\fixedpointop$, \[ \weights^{(\iteridx+1)} = \fixedpointop\big(\weights^{(\iteridx)}\big) \text{, for } \iteridx = 0,1,\ldots \text{.} \] The operator $\fixedpointop$ is designed such that the sequence of model parameters $\weights^{(\iteridx)}$ converges to a solution of the normal equations \eqref{eq_linreg_normal_eq_dict}. \par One such iterative optimization method is gradient descent (GD), which uses the GD step operator \[ \gdstep{\lrate}(\weights) = \weights - \lrate \nabla f(\weights) \text{.} \] Here, $\lrate > 0$ is a step size that must be chosen small enough to ensure convergence (Bertsekas, 2016, Prop.~1.2.2). The superscript in $\gdstep{\lrate}$ makes the dependence of the operator on $\lrate$ explicit. Inserting the average training error into the general form of the GD step yields the explicit update \[ \begin{aligned} \gdstep{\lrate}(\weights) &= \weights + \frac{2\lrate}{\samplesize} \sum_{\sampleidx=1}^{\samplesize} \big( \truelabel^{(\sampleidx)} - \weights^{\top}\featurevec^{(\sampleidx)} \big) \featurevec^{(\sampleidx)} \\ &= \Big(\mI - \frac{2\lrate}{\samplesize}\featuremtx^{\top}\featuremtx\Big)\weights + \frac{2\lrate}{\samplesize}\featuremtx^{\top}\labelvec \text{.} \end{aligned} \] \par The fixed points of $\gdstep{\lrate}$ are exactly the solutions of the normal equations \eqref{eq_linreg_normal_eq_dict}, since $\gdstep{\lrate}(\widehat{\weights}) = \widehat{\weights}$ holds if and only if $\featuremtx^{\top}\featuremtx\widehat{\weights} = \featuremtx^{\top}\labelvec$. The update operator is affine, $$\gdstep{\lrate}(\weights) = \mM \weights + \vb \text{,}$$ with linear part $\mM = \mI - (2\lrate/\samplesize)\featuremtx^{\top}\featuremtx$ and offset $\vb = (2\lrate/\samplesize)\, \featuremtx^{\top}\labelvec$. \par When $\featuremtx$ has full column rank, so that $\featuremtx^{\top}\featuremtx$ is positive definite (pd), $\gdstep{\lrate}$ is a contractive operator with respect to the Euclidean norm for every step size $0 < \lrate < \samplesize/\eigval{\max}$. Here, $\eigval{\max}$ is the largest eigenvalue of $\featuremtx^{\top}\featuremtx$. The eigenvalues of the linear part $\mM$ are $1 - (2\lrate/\samplesize)\,\eigval{\featureidx}$, with the positive eigenvalues $\eigval{\featureidx}$ of $\featuremtx^{\top}\featuremtx$. They lie strictly between $-1$ and $1$ exactly for such step sizes. The more general case of any smooth, strongly convex objective function and step sizes $\lrate \leq \samplesize/(\eigval{\min} + \eigval{\max})$ is covered by Nesterov (2004), Thm.~2.1.15, with the smallest eigenvalue $\eigval{\min}$ of $\featuremtx^{\top}\featuremtx$. The fixed-point iteration converges to the unique solution $\widehat{\weights}$. Moreover, the convergence speed of GD is governed by the ratio of the largest to the smallest eigenvalue of $\featuremtx^{\top}\featuremtx$, that is, by its condition number $\condnumber{\featuremtx^{\top}\featuremtx} = \eigval{\max}/\eigval{\min}$ (Boyd and Vandenberghe, 2004, Sect.~9.3). \par The GD step above uses the entire training set to compute the gradient. When the data points instead arrive sequentially at time instants $\iteridx=1,2,\ldots$, or the training set is too large to fit in memory, the model parameters can be learned by online gradient descent (online GD). At each time step $\iteridx$, it applies a single GD step. That step uses only the arriving data point $\pair{\featurevec^{(\iteridx)}}{\truelabel^{(\iteridx)}}$, \[ \weights^{(\iteridx+1)} = \weights^{(\iteridx)} + 2\lrate \big( \truelabel^{(\iteridx)} - \big(\weights^{(\iteridx)}\big)^{\top} \featurevec^{(\iteridx)} \big) \featurevec^{(\iteridx)} \text{.} \] This is the GD step on the loss function $\big(\truelabel^{(\iteridx)} - \weights^{\top}\featurevec^{(\iteridx)}\big)^{2}$ of the arriving data point alone, rather than on the average squared error loss $f(\weights)$. It is the least mean squares (LMS) update (Bishop, 2006, Sect.~3.1.3). It avoids storing the entire training set, and its per-iteration complexity is independent of the training set size $\samplesize$. \par The optimality condition \eqref{eq_linreg_normal_eq_dict} allows one to study the stability of linear regression. Ideally, the learned model parameters are insensitive to limited perturbations of the training set. A perturbed training set yields a feature matrix $\widetilde{\featuremtx} = \featuremtx + \Delta \featuremtx$ and a label vector $\widetilde{\labelvec} = \labelvec + \Delta \labelvec$, with perturbation matrix $\Delta \featuremtx$ and vector $\Delta \labelvec$. This gives the perturbed normal equations \begin{equation} \label{eq_linreg_perturbed_normal_eq_dict} \widetilde{\featuremtx}^{\top} \widetilde{\featuremtx} \widetilde{\weights} = \widetilde{\featuremtx}^{\top} \widetilde{\labelvec} \text{.} \end{equation} Matrix perturbation theory quantifies how much $\widetilde{\weights}$ deviates from a solution $\widehat{\weights}$ of the clean normal equations \eqref{eq_linreg_normal_eq_dict} (Golub and Van Loan, 2013, Sect.~5.3.6). \par For a label-only perturbation, which replaces a single label of the training set with an outlier or another corrupted value, $\Delta \featuremtx = \mathbf{0}$, and the perturbed and clean normal equations share the coefficient matrix $\featuremtx^{\top}\featuremtx$. Assume $\featuremtx$ has full column rank, so $\featuremtx^{\top}\featuremtx$ is invertible and the solution is unique. Subtracting the clean normal equations \eqref{eq_linreg_normal_eq_dict} from the perturbed normal equations \eqref{eq_linreg_perturbed_normal_eq_dict} yields \[ \widetilde{\weights} - \widehat{\weights} = \big(\featuremtx^{\top}\featuremtx\big)^{-1} \featuremtx^{\top} \Delta \labelvec = \pinv{\featuremtx} \Delta \labelvec \text{,} \] with the pseudoinverse $\pinv{\featuremtx}$. The effect of perturbing the training set is thereby quantified, \[ \normgeneric{\widetilde{\weights} - \widehat{\weights}}{2} \leq \normgeneric{\pinv{\featuremtx}}{2}\, \normgeneric{\Delta \labelvec}{2} \text{.} \] When $\featuremtx$ has full column rank, the pseudoinverse has the norm $\normgeneric{\pinv{\featuremtx}}{2} = 1/\sqrt{\eigval{\min}(\featuremtx^{\top}\featuremtx)} = \sqrt{\condnumber{\featuremtx^{\top}\featuremtx}}/\normgeneric{\featuremtx}{2}$, where $\condnumber{\featuremtx^{\top}\featuremtx} = \eigval{\max}/\eigval{\min}$ is the condition number of $\featuremtx^{\top}\featuremtx$ (see pseudoinverse). In the setting of Fig.~\ref{fig_linreg_dict}, $\featuremtx = \mathbf{1}$ (the all-ones vector), so $\pinv{\featuremtx} \Delta \labelvec = (1/\samplesize) \sum_{\sampleidx=1}^{\samplesize} \Delta \truelabel^{(\sampleidx)}$. Fig.~\ref{fig_linreg_outlier_dict} illustrates such a label-only perturbation of the training set from Fig.~\ref{fig_linreg_dict}: a single label perturbation shifts $\widehat{w}$ by one third of it, as Fig.~\ref{fig_linreg_outlier_dict} shows. \begin{figure}[H] \centering \begin{tikzpicture}[scale=0.7, y=1.1cm, x=0.5cm] \foreach \x/\y in { 1/2, 4/3 } { \filldraw[blue] (\x, \y) circle (2pt); } \filldraw[blue] (7, 4) circle (2pt); \draw[blue, thick] (7, 10) circle (2pt); \draw[->, red, thick] (7, 4.5) -- (7, 9.5) node[midway, right] {$\Delta \truelabel^{(3)} = 6$}; \node[above left=2pt and 2pt, red] at (7, 10) {outlier}; \draw[dashed, thick] (0.5, 3) -- (10.5, 3) node[right] {$\widehat{w}=3$}; \draw[dotted, thick] (0.5, 5) -- (10.5, 5) node[right] {$\widetilde{w}=5$}; \end{tikzpicture} \caption{The perturbation $\Delta \truelabel^{(3)} = 6$ (red arrow) replaces the third label of the training set from Fig.~\ref{fig_linreg_dict} with an outlier: filled circles denote the original data points, as in Fig.~\ref{fig_linreg_dict}, and the open circle the perturbed value. The outlier raises the average of the labels: the learned model parameter shifts by $\widetilde{w} - \widehat{w} = 6/3 = 2$, from the clean solution $\widehat{w} = 3$ (dashed line) to the perturbed solution $\widetilde{w} = 5$ (dotted line) \label{fig_linreg_outlier_dict}} \end{figure} The effect of perturbations can also be studied on the level of a specific optimization method. For example, the update of online GD becomes \[ \begin{aligned} \weights^{(\iteridx+1)} &= \weights^{(\iteridx)} + 2\lrate \big( \widetilde{\truelabel}^{(\iteridx)} - \big(\weights^{(\iteridx)}\big)^{\top} \widetilde{\featurevec}^{(\iteridx)} \big) \widetilde{\featurevec}^{(\iteridx)} \\ &= \weights^{(\iteridx)} + 2\lrate \big( \truelabel^{(\iteridx)} - \big(\weights^{(\iteridx)}\big)^{\top} \featurevec^{(\iteridx)} \big) \featurevec^{(\iteridx)} + \perturbation{\iteridx} \text{.} \end{aligned} \] Here, $\perturbation{\iteridx}$ is a perturbation term that depends on the data point perturbation $\pair{\Delta \featurevec^{(\iteridx)}}{\Delta \truelabel^{(\iteridx)}}$ and the current model parameters $\weights^{(\iteridx)}$. A bounded perturbation term moves each iterate by a bounded amount. The perturbed online GD update has no fixed point in general, since a single perturbed data point leaves $\widetilde{\weights}$ unchanged only if $\widetilde{\weights}$ interpolates it. The GD step operator $\gdstep{\lrate}$ computed on the perturbed training set instead has as its fixed points exactly the solutions of the perturbed normal equations \eqref{eq_linreg_perturbed_normal_eq_dict}, by the same computation as for the clean training set. Among them is the perturbed solution $\widetilde{\weights}$, whose displacement from $\widehat{\weights}$ was computed above. \subsection{Cross References} regression, linear model, least squares, squared error loss, ERM, GD, ridge regression, Lasso, pseudoinverse \subsection{References} \begin{enumerate} \item Bertsekas (2016). Nonlinear Programming. 3rd ed., Athena Scientific, Belmont, MA, USA. \item Bishop (2006). Pattern Recognition and Machine Learning. Springer Science+Business Media, New York, NY, USA. \item Boyd and Vandenberghe (2004). Convex Optimization. Cambridge Univ. Press, Cambridge, U.K. \item Golub and Van Loan (2013). Matrix Computations. 4th ed., The Johns Hopkins Univ. Press, Baltimore, MD, USA. \item Pedregosa et al. (2011). Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research, 12(85), pp.~2825--2830. \item Lehmann and Casella (1998). Theory of Point Estimation. 2nd ed., Springer-Verlag, New York, NY, USA. \item Meyer (2000). Matrix Analysis and Applied Linear Algebra. SIAM, Philadelphia, PA, USA. \item Paszke et al. (2019). PyTorch: An Imperative Style, High-Performance Deep Learning Library. In: Adv. Neural Inf. Process. Syst., pp.~8026--8037. Curran Associates Inc. \item Strang (2016). Introduction to Linear Algebra. 5th ed., Wellesley-Cambridge Press, Wellesley, MA, USA. \item Hastie et al. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2nd ed., Springer Science+Business Media, New York, NY, USA. \item Nesterov (2004). Introductory Lectures on Convex Optimization: A Basic Course. Kluwer Academic, Boston, MA, USA. \item Papoulis and Pillai (2002). Probability, Random Variables, and Stochastic Processes. 4th ed., McGraw-Hill Higher Education, New York, NY, USA. \end{enumerate} \end{document}