Linear models¶

Basics of modeling, optimization, and regularization

Mahmood Amintoosi, Spring 2026

Computer Science Dept, Ferdowsi University of Mashhad

I should mention that the original material of this course was from Open Machine Learning Course, by Joaquin Vanschoren and others.

Notation and Definitions¶

• A scalar is a simple numeric value, denoted by an italic letter: $x=3.24$

• A vector is a 1D ordered array of n scalars, denoted by a bold letter: $\mathbf{x}=[3.24, 1.2]$

  • $x_i$ denotes the $i$th element of a vector, thus $x_0 = 3.24$.

    • Note: some other courses use $x^{(i)}$ notation

• A set is an unordered collection of unique elements, denote by caligraphic capital: $\mathcal{S}=\{3.24, 1.2\}$

• A matrix is a 2D array of scalars, denoted by bold capital: $\mathbf{X}=\begin{bmatrix} 3.24 & 1.2 \\ 2.24 & 0.2 \end{bmatrix}$

  • $\textbf{X}_{i}$ denotes the $i$th row of the matrix

  • $\textbf{X}_{:,j}$ denotes the $j$th column

  • $\textbf{X}_{i,j}$ denotes the element in the $i$th row, $j$th column, thus $\mathbf{X}_{1,0} = 2.24$

• $\mathbf{X}^{n \times p}$, an $n \times p$ matrix, can represent $n$ data points in a $p$-dimensional space

  • Every row is a vector that can represent a point in an p-dimensional space, given a basis.

  • The standard basis for a Euclidean space is the set of unit vectors

• E.g. if $\mathbf{X}=\begin{bmatrix} 3.24 & 1.2 \\ 2.24 & 0.2 \\ 3.0 & 0.6 \end{bmatrix}$

• A tensor is an k-dimensional array of data, denoted by an italic capital: $T$

  • k is also called the order, degree, or rank

  • $T_{i,j,k,...}$ denotes the element or sub-tensor in the corresponding position

  • A set of color images can be represented by:

    • a 4D tensor (sample x height x width x color channel)

    • a 2D tensor (sample x flattened vector of pixel values)

Basic operations¶

• Sums and products are denoted by capital Sigma and capital Pi:

• Operations on vectors are element-wise: e.g. $\mathbf{x}+\mathbf{z} = [x_0+z_0,x_1+z_1, ... , x_p+z_p]$

• Dot product $\mathbf{w}\mathbf{x} = \mathbf{w} \cdot \mathbf{x} = \mathbf{w}^{T} \mathbf{x} = \sum_{i=0}^{p} w_i \cdot x_i = w_0 \cdot x_0 + w_1 \cdot x_1 + ... + w_p \cdot x_p$

• Matrix product $\mathbf{W}\mathbf{x} = \begin{bmatrix} \mathbf{w_0} \cdot \mathbf{x} \\ ... \\ \mathbf{w_p} \cdot \mathbf{x} \end{bmatrix}$

• A function $f(x) = y$ relates an input element $x$ to an output $y$

  • It has a local minimum at $x=c$ if $f(x) \geq f(c)$ in interval $(c-\epsilon, c+\epsilon)$

  • It has a global minimum at $x=c$ if $f(x) \geq f(c)$ for any value for $x$

• A vector function consumes an input and produces a vector: $\mathbf{f}(\mathbf{x}) = \mathbf{y}$

• $\underset{x\in X}{\operatorname{max}}f(x)$ returns the largest value f(x) for any x

• $\underset{x\in X}{\operatorname{argmax}}f(x)$ returns the element x that maximizes f(x)

Gradients¶

• A derivative $f'$ of a function $f$ describes how fast $f$ grows or decreases

• The process of finding a derivative is called differentiation

  • Derivatives for basic functions are known

  • For non-basic functions we use the chain rule: $F(x) = f(g(x)) \rightarrow F'(x)=f'(g(x))g'(x)$

• A function is differentiable if it has a derivative in any point of it’s domain

  • It’s continuously differentiable if $f'$ is a continuous function

  • We say $f$ is smooth if it is infinitely differentiable, i.e., $f', f'', f''', ...$ all exist

• A gradient $\nabla f$ is the derivative of a function in multiple dimensions

  • It is a vector of partial derivatives: $\nabla f = \left[ \frac{\partial f}{\partial x_0}, \frac{\partial f}{\partial x_1},... \right]$

  • E.g. $f=2x_0+3x_1^{2}-\sin(x_2) \rightarrow \nabla f= [2, 6x_1, -cos(x_2)]$

• Example: $f = -(x_0^2+x_1^2)$

  • $\nabla f = \left[\frac{\partial f}{\partial x_0},\frac{\partial f}{\partial x_1}\right] = \left[-2x_0,-2x_1\right]$

  • Evaluated at point (-4,1): $\nabla f(-4,1) = [8,-2]$

    • These are the slopes at point (-4,1) in the direction of $x_0$ and $x_1$ respectively

Distributions and Probabilities¶

• The normal (Gaussian) distribution with mean $\mu$ and standard deviation $\sigma$ is noted as $N(\mu,\sigma)$

• A random variable $X$ can be continuous or discrete

• A probability distribution $f_X$ of a continuous variable $X$: probability density function (pdf)

  • The expectation is given by $\mathbb{E}[X] = \int x f_{X}(x) dx$

• A probability distribution of a discrete variable: probability mass function (pmf)

  • The expectation (or mean) $\mu_X = \mathbb{E}[X] = \sum_{i=1}^k[x_i \cdot Pr(X=x_i)]$

Linear models¶

Linear models make a prediction using a linear function of the input features $X$

Learn $w$ from $X$, given a loss function $\mathcal{L}$:

• Many algorithms with different $\mathcal{L}$: Least squares, Ridge, Lasso, Logistic Regression, Linear SVMs,...

• Can be very powerful (and fast), especially for large datasets with many features.

• Can be generalized to learn non-linear patterns: Generalized Linear Models

  • Features can be augmentented with polynomials of the original features

  • Features can be transformed according to a distribution (Poisson, Tweedie, Gamma,...)

  • Some linear models (e.g. SVMs) can be kernelized to learn non-linear functions

Linear models for regression¶

• Prediction formula for input features x:

  • $w_1$ ... $w_p$ usually called weights or coefficients , $w_0$ the bias or intercept

  • Assumes that errors are $N(0,\sigma)$

Linear Regression (aka Ordinary Least Squares)¶

• Loss function is the sum of squared errors (SSE) (or residuals) between predictions $\hat{y}_i$ (red) and the true regression targets $y_i$ (blue) on the training set.

Solving ordinary least squares¶

Closed-Form Solution

• Convex optimization problem with unique closed-form solution:

  $$\mathbf{w}^{*} = (X^{T}X)^{-1} X^T \mathbf{y}$$

  (6)

If we have:

then $\hat{y}$ can be considered as dot product of two vectors:

If $i^{th}$ instance is shown by $\mathbf{x}_i$

then:

Loss: Sum of Squared Error:

where:

We Know that (See Matrix Calculus in Wikipedia):

and if $M$ be a symmetric matrix:

Hence:

Set the derivative equal to zero:

Some notes about Closed-Form Solution¶

The unique closed-form solution:

Where X has $n$ rows, $p$ features, hence $X^{T}X$ has dimensionality $p \cdot p$

Notes:

• Slow. Time complexity is quadratic in number of features: $\mathcal{O}(p^2n)$

• Only works if $n>p$, $X^T X$ is Singular When $n < p$, making OLS’s $(X^T X)^{-1}$ undefined.

• The closed-form solution of Ordinary Least Squares (OLS) is highly prone to overfitting, especially when applied to small or highly correlated datasets.

• Large coefficients: In the absence of regularization, weights $w$ can grow excessively, leading to unstable predictions.

• Sensitive to input variations: Small changes in $x$ may result in disproportionately large shifts in the output $y$, making the model unreliable.

• No direct hyperparameter control: Unlike Gradient Descent, OLS lacks tuning mechanisms such as learning rate or regularization strength.

Gradient Descent¶

• More efficient for large-scale or high-dimensional datasets.

• Preferred when computing $X^{T}X$ is infeasible or excessively time-consuming due to large $p$ or $n$.

• Offers multiple tunable settings to control the learning process, including:

  • Learning Rate: Adjusts update speed to balance stability and convergence.

  • Regularization: Options like L1, L2, or Elastic Net mitigate overfitting.

  • Batch Size: Influences computational efficiency (Stochastic, Mini-batch, Full-batch).

  • Iterations: Determines the number of passes over the dataset.

  • Momentum & Decay: Helps stabilize and refine the learning trajectory.

Gradient Descent¶

• Start with an initial, random set of weights: $\mathbf{w}^0$

• Given a differentiable loss function $\mathcal{L}$ (e.g. $\mathcal{L}_{SSE}$), compute $\nabla \mathcal{L}$

• For least squares: $\frac{\partial \mathcal{L}_{SSE}}{\partial w_i}(\mathbf{w}) = -2\sum_{n=1}^{N} (y_n-\hat{y}_n) x_{n,i}$

  • If feature $X_{:,i}$ is associated with big errors, the gradient wrt $w_i$ will be large

• Update all weights slightly (by step size or learning rate $\eta$) in ‘downhill’ direction.

• Basic update rule (step s):

  $$\mathbf{w}^{s+1} = \mathbf{w}^s-\eta\nabla \mathcal{L}(\mathbf{w}^s)$$

  (20)

• Important hyperparameters

  • Learning rate

    • Too small: slow convergence. Too large: possible divergence

  • Maximum number of iterations

    • Too small: no convergence. Too large: wastes resources

  • Learning rate decay with decay rate $k$

    • E.g. exponential ($\eta^{s+1} = \eta^{0} e^{-ks}$), inverse-time ($\eta^{s+1} = \frac{\eta^{s}}{1+ks}$),...

  • Many more advanced ways to control learning rate (see later)

    • Adaptive techniques: depend on how much loss improved in previous step

Effect of learning rate

Effect of learning rate decay

In two dimensions:

• You can get stuck in local minima (if the loss is not fully convex)

  • If you have many model parameters, this is less likely

  • You always find a way down in some direction

  • Models with many parameters typically find good local minima

• Intuition: walking downhill using only the slope you “feel” nearby

(Image by A. Karpathy)

Stochastic Gradient Descent (SGD)¶

• Compute gradients not on the entire dataset, but on a single data point $i$ at a time

  • Gradient descent: $\mathbf{w}^{s+1} = \mathbf{w}^s-\eta\nabla \mathcal{L}(\mathbf{w}^s) = \mathbf{w}^s-\frac{\eta}{n} \sum_{i=1}^{n} \nabla \mathcal{L_i}(\mathbf{w}^s)$

  • Stochastic Gradient Descent: $\mathbf{w}^{s+1} = \mathbf{w}^s-\eta\nabla \mathcal{L_i}(\mathbf{w}^s)$

• Many smoother variants, e.g.

  • Minibatch SGD: compute gradient on batches of data: $\mathbf{w}^{s+1} = \mathbf{w}^s-\frac{\eta}{B} \sum_{i=1}^{B} \nabla \mathcal{L_i}(\mathbf{w}^s)$

  • Stochastic Average Gradient Descent (SAG). With $i_s \in [1,n]$ randomly chosen per iteration:

    • Incremental gradient: $\mathbf{w}^{s+1} = \mathbf{w}^s-\frac{\eta}{n} \sum_{i=1}^{n} v_i^s$ with $v_i^s = \begin{cases}\nabla \mathcal{L_i}(\mathbf{w}^s) & i = i_s \\ v_i^{s-1} & \text{otherwise} \end{cases}$

In practice¶

• Linear regression can be found in sklearn.linear_model. We’ll evaluate it on the Boston Housing dataset.

  • LinearRegression uses closed form solution, SGDRegressor with loss='squared_loss' uses Stochastic Gradient Descent

  • Large coefficients signal overfitting

  • Test score is much lower than training score

Ridge regression¶

• Adds a penalty term to the least squares loss function:

• Model is penalized if it uses large coefficients ($w$)

  • Each feature should have as little effect on the outcome as possible

  • We don’t want to penalize $w_0$, so we leave it out

• Regularization: explicitly restrict a model to avoid overfitting.

  • Called L2 regularization because it uses the L2 norm: $\sum w_i^2$

• The strength of the regularization can be controlled with the $\alpha$ hyperparameter.

  • Increasing $\alpha$ causes more regularization (or shrinkage). Default is 1.0.

• Still convex. Can be optimized in different ways:

  • Closed form solution (a.k.a. Cholesky): $w^{*} = (X^{T}X + \alpha I)^{-1} X^T \mathbf{y}$

  • Gradient descent and variants, e.g. Stochastic Average Gradient (SAG,SAGA)

    • Conjugate gradient (CG): each new gradient is influenced by previous ones

  • Use Cholesky for smaller datasets, Gradient descent for larger ones

Ridge Regression Derivation¶

We want to minimize the Ridge loss function with respect to $\mathbf{w}$:

Expanding:

Taking the derivative with respect to $\mathbf{w}$:

Setting the derivative equal to zero:

In practice¶

from sklearn.linear_model import Ridge
lr = Ridge().fit(X_train, y_train)

Test set score is higher and training set score lower: less overfitting!

• We can plot the weight values for differents levels of regularization to explore the effect of $\alpha$.

• Increasing regularization decreases the values of the coefficients, but never to 0.

• When we plot the train and test scores for every $\alpha$ value, we see a sweet spot around $\alpha=0.2$

  • Models with smaller $\alpha$ are overfitting

  • Models with larger $\alpha$ are underfitting

Other ways to reduce overfitting¶

• Add more training data: with enough training data, regularization becomes less important

  • Ridge and ordinary least squares will have the same performance

• Use fewer features: remove unimportant ones or find a low-dimensional embedding (e.g. PCA)

  • Fewer coefficients to learn, reduces the flexibility of the model

• Scaling the data typically helps (and changes the optimal $\alpha$ value)

Lasso (Least Absolute Shrinkage and Selection Operator)¶

• Adds a different penalty term to the least squares sum:

• Called L1 regularization because it uses the L1 norm

  • Will cause many weights to be exactly 0

• Same parameter $\alpha$ to control the strength of regularization.

  • Will again have a ‘sweet spot’ depending on the data

• No closed-form solution

• Convex, but no longer strictly convex, and not differentiable

  • Weights can be optimized using coordinate descent

Analyze what happens to the weights:

• L1 prefers coefficients to be exactly zero (sparse models)

• Some features are ignored entirely: automatic feature selection

• How can we explain this?

Coordinate descent¶

• Alternative for gradient descent, supports non-differentiable convex loss functions (e.g. $\mathcal{L}_{Lasso}$)

• In every iteration, optimize a single coordinate $w_i$ (find minimum in direction of $x_i$)

  • Continue with another coordinate, using a selection rule (e.g. round robin)

• Faster iterations. No need to choose a step size (learning rate).

• May converge more slowly. Can’t be parallellized.

Coordinate descent with Lasso (Further Reading)¶

• Remember that $\mathcal{L}_{Lasso} = \mathcal{L}_{SSE} + \alpha \sum_{i=1}^{p} |w_i|$

• For one $w_i$: $\mathcal{L}_{Lasso}(w_i) = \mathcal{L}_{SSE}(w_i) + \alpha |w_i|$

• The L1 term is not differentiable but convex: we can compute the subgradient

  • Unique at points where $\mathcal{L}$ is differentiable, a range of all possible slopes [a,b] where it is not

  • For $|w_i|$, the subgradient $\partial_{w_i} |w_i|$ = $\begin{cases}-1 & w_i<0\\ [-1,1] & w_i=0 \\ 1 & w_i>0 \\ \end{cases}$

  • Subdifferential $\partial(f+g) = \partial f + \partial g$ if $f$ and $g$ are both convex

• To find the optimum for Lasso $w_i^{*}$, solve

  $$\begin{aligned} \partial_{w_i} \mathcal{L}_{Lasso}(w_i) &= \partial_{w_i} \mathcal{L}_{SSE}(w_i) + \partial_{w_i} \alpha |w_i| \\ 0 &= (w_i - \rho_i) + \alpha \cdot \partial_{w_i} |w_i| \\ w_i &= \rho_i - \alpha \cdot \partial_{w_i} |w_i| \end{aligned}$$

  (27)

  • In which $\rho_i$ is the part of $\partial_{w_i} \mathcal{L}_{SSE}(w_i)$ excluding $w_i$ (assume $z_i=1$ for now)

    • $\rho_i$ can be seen as the $\mathcal{L}_{SSE}$ ‘solution’: $w_i = \rho_i$ if $\partial_{w_i} \mathcal{L}_{SSE}(w_i) = 0$

  $$\partial_{w_i} \mathcal{L}_{SSE}(w_i) = \partial_{w_i} \sum_{n=1}^{N} (y_n-(\mathbf{w}\mathbf{x_n} + w_0))^2 = z_i w_i -\rho_i$$

  (28)

• We found: $w_i = \rho_i - \alpha \cdot \partial_{w_i} |w_i|$

• The Lasso solution has the form of a soft thresholding function $S$

  $$w_i^* = S(\rho_i,\alpha) = \begin{cases} \rho_i + \alpha, & \rho_i < -\alpha \\ 0, & -\alpha < \rho_i < \alpha \\ \rho_i - \alpha, & \rho_i > \alpha \\ \end{cases}$$

  (29)

  • Small weights become 0: sparseness!

  • If the data is not normalized, $w_i^* = \frac{1}{z_i}S(\rho_i,\alpha)$ with constant $z_i = \sum_{n=1}^{N} x_{ni}^2$

• Ridge solution: $w_i = \rho_i - \alpha \cdot \partial_{w_i} w_i^2 = \rho_i - 2\alpha \cdot w_i$, thus $w_i^* = \frac{\rho_i}{1 + 2\alpha}$

• For mor information see Rashidabadi (2016)

Interpreting L1 and L2 loss¶

• In 2D (for 2 model weights $w_1$ and $w_2$)

  • The least squared loss is a 2D convex function in this space

  • For illustration, assume that L1 loss = L2 loss = 1

    • L1 loss ($\Sigma |w_i|$): every {$w_1, w_2$} falls on the diamond

    • L2 loss ($\Sigma w_i^2$): every {$w_1, w_2$} falls on the circle

  • For L1, the loss is minimized if $w_1$ or $w_2$ is 0 (rarely so for L2)

Comparing Ridge and LASSO for Prediction and Feature Selection¶

The following code trains and compares Linear Regression, Ridge (L2), and LASSO (L1) regression models. We’ll examine their predictive performance using Mean Squared Error (MSE) and observe how LASSO can be used for feature selection by shrinking some feature coefficients to zero.

Elastic-Net¶

• Adds both L1 and L2 regularization:

• $\rho$ is the L1 ratio

  • With $\rho=1$, $\mathcal{L}_{Elastic} = \mathcal{L}_{Lasso}$

  • With $\rho=0$, $\mathcal{L}_{Elastic} = \mathcal{L}_{Ridge}$

  • $0 < \rho < 1$ sets a trade-off between L1 and L2.

• Allows learning sparse models (like Lasso) while maintaining L2 regularization benefits

  • E.g. if 2 features are correlated, Lasso likely picks one randomly, Elastic-Net keeps both

• Weights can be optimized using coordinate descent (similar to Lasso)

Sparse Optimization: Theory and Applications¶

Sparse optimization is a fundamental paradigm in mathematical modeling that seeks solutions with few non-zero elements, often formalized through ℓ₁-norm regularization or non-convex penalties. This framework is grounded in the Compressed Sensing (CS) theory Donoho (2006), which guarantees exact signal recovery from sub-Nyquist measurements if the signal is sparse in some basis. Key applications include:

1. Medical Imaging: CS enables faster MRI scans by reconstructing images from limited k-space samples (Lustig et al., 2007), as highlighted in the CUHK lecture notes on CS (p. 13).

2. Signal Processing: Sparse methods underpin denoising (e.g., wavelet shrinkage) and inpainting, where total variation minimization preserves edges while removing noise (Figueiredo, 2021, slides p. 8).

3. Machine Learning: Sparse regression (e.g., LASSO; Tibshirani, 1996) and robust PCA (Candès et al., 2011) are critical for feature selection and anomaly detection. The Ghent tutorial (p. 5) further links sparsity to deep network pruning.

4. Statistics: High-dimensional inference (e.g., genomics) benefits from sparsity-induced interpretability (Bühlmann & Van de Geer, 2011).

The field continues to evolve with non-convex penalties (e.g., SCAD) and greedy algorithms (OMP), balancing computational efficiency and statistical guarantees. For a unified perspective, see the cited references and tutorials.

Sparse Signal Denoising via Weighted Lasso in DCT Domain¶

Given a noisy 1D signal:

where:

• f(t) is the clean signal (sparse in frequency domain)

• ϵ(t) is Gaussian noise (non-sparse in any basis)

We solve the weighted Lasso problem:

Noise-Sparsity Dichotomy¶

• Signal: Sparse in DCT domain (few dominant coefficients)

• Noise: Non-sparse (energy spread across all frequencies)

Weighted ℓ₁ Magic¶

• Preserves signal structure (low frequencies)

• Aggressively truncates noise (high frequencies)

Performance Analysis¶

Theoretical Guarantees¶

• Exact Recovery: Possible if:

  • Signal truly sparse in DCT basis

  • Noise level below threshold (Donoho et al., 2006)

• Stability: Weighted Lasso provides better coefficient preservation than standard Lasso

Limitations¶

• Assumes signal sparsity in chosen basis

• Requires tuning for each application

• Computationally heavier than simple filtering

Signal Denoising using LASSO

Orthogonal Matching Pursuit (OMP) for Sparse Recovery¶

Orthogonal Matching Pursuit (OMP) is a greedy algorithm that solves the sparse approximation problem:

Problem Statement: Given a design matrix $X \in \mathbb{R}^{n \times p}$ (where typically $n \ll p$) and observations $\mathbf{y} \in \mathbb{R}^n$, find the sparsest vector $\mathbf{w} \in \mathbb{R}^p$ such that $\mathbf{y} \approx X\mathbf{w}$.

Key characteristics:

• Iterative approach that selects one feature per iteration based on maximum correlation

• Provides exact recovery guarantees under Restricted Isometry Property (RIP) conditions

• Computationally efficient ($O(knp)$ complexity) compared to $\ell_0$-norm minimization

• Widely used in compressed sensing, signal processing, and feature selection

• Preserves interpretability through explicit feature selection

This method complements LASSO by offering an alternative approach to sparse recovery with distinct theoretical guarantees and computational trade-offs.

While LASSO solves the $\ell_1$-regularized problem:

OMP provides a greedy alternative for sparse approximation, solving:

OMP algorithm¶

Given:

• Design matrix $X \in \mathbb{R}^{n \times p}$

• Response vector $\mathbf{y} \in \mathbb{R}^n$

• Target sparsity level $k$

At each iteration:

1. Select most correlated feature: $j = \arg\max_{j \notin \mathcal{A}} |\mathbf{x}_j^\top \mathbf{r}|$, where $\mathbf{r}$ is the current residual and $\mathcal{A}$ is the active feature set.

2. Update support: $\mathcal{A} \leftarrow \mathcal{A} \cup \{j\}$

3. Solve least squares: $\mathbf{w} = \arg\min_{\mathbf{w}} \|X_{\mathcal{A}}\mathbf{w} - \mathbf{y}\|_2^2$

4. Update residual: $\mathbf{r} \leftarrow X_{\mathcal{A}}\mathbf{w} - \mathbf{y}$

Termination: After $k$ iterations or when $\|\mathbf{r}\|_2 < \epsilon$

Figure 10.2 of Theodoridis (2020):

(A) In the case of an orthogonal matrix, the observations vector y will be orthogonal to any inactive column, here, $x_3^c$. (B) In the more general case, it is expected to “lean” closer (form smaller angles) to the active than to the inactive columns.

Figure 2.3 of Bakhshali (2018)

Comparison with LASSO¶

For further information see Asaran (2016).

Practice¶

A Simple Example

When to Use OMP¶

1. Exactly sparse signals – OMP is ideal when the underlying signal truly has only a few non-zero coefficients, making it a strong choice for enforcing strict sparsity.

2. Small-scale problems – Works efficiently when the number of non-zero elements $k$ is much smaller than the total features $p$ (i.e., $k \ll p$), preventing unnecessary computational overhead.

3. Interpretable solutions – Since OMP selects features iteratively, the stepwise approach makes it easier to analyze which features contribute to the model, often simpler than interpreting LASSO’s continuous shrinkage.

4. Dictionary learning – Frequently used in sparse coding applications, where signals are represented as a sparse combination of dictionary elements, making it valuable for feature selection in compressed sensing.

Similarities and Differences in Solving $A\mathbf{x} = \mathbf{b}$ in Regression and Sparse Representation¶

Similarities¶

1. Core Equation: Both domains solve the linear system $A\mathbf{x} = \mathbf{b}$.

2. Matrix-Vector Structure:

   • $A$ is a matrix, $\mathbf{x}$ is a vector of unknowns, and $\mathbf{b}$ is a vector of observations/targets.

3. Optimization Goal: Both aim to find $\mathbf{x}$ such that $A\mathbf{x}$ approximates $\mathbf{b}$, often with a minimization objective (e.g., least squares or sparsity constraints).

Differences in Specification of $A, \mathbf{x}, \mathbf{b}$¶

Key Differences Explained¶

1. Structure of $A$:

   • Regression: $A$ is tall/square (rows ≥ columns). Rows = samples, columns = features.
     Example: 100 patients (rows) × 5 features (columns).

   • Sparse Representation: $A$ is fat (columns > rows). Columns = atoms, rows = signal dimensions.
     Example: 1000 pixels (rows) × 5000 atoms (columns).

2. Sparsity of $\mathbf{x}$:

   • Regression: $\mathbf{x}$ is dense (all entries non-zero). Every feature contributes to predictions.

   • Sparse Representation: $\mathbf{x}$ must be sparse (most entries zero). Only a few atoms are activated.

3. Interpretation:

   • Regression: $\mathbf{x}$ quantifies feature importance (e.g., “size” has weight 0.8).

   • Sparse Representation: $\mathbf{x}$ identifies active atoms (e.g., “edge #42” is activated for a line in an image).

4. Optimization:

   • Regression: Unconstrained least squares (closed-form or iterative solvers).

   • Sparse Representation: Constrained optimization (L1/L0 regularization) to enforce sparsity (e.g., LASSO, OMP).

Example: House Price Prediction (Regression)¶

• $A$: 100 houses (rows) × 3 features: [size, bedrooms, age].

• $\mathbf{b}$: House prices (e.g., [300k, 450k, ...]).

• $\mathbf{x}$: Weights (e.g., [size: 0.2, bedrooms: 0.1, age: -0.05]).

• Solution: $\mathbf{x} = (A^T A)^{-1} A^T \mathbf{b}$ (dense).

Example: Image Denoising (Sparse Representation)¶

• $A$: Dictionary of 10,000 image patches (atoms).

• $\mathbf{b}$: Noisy 50×50 image vectorized (2500 rows).

• $\mathbf{x}$: Sparse vector activating 20 atoms (e.g., only 20 non-zero entries).

• Solution: Minimize $\|A\mathbf{x} - \mathbf{b}\|_2^2 + \lambda \|\mathbf{x}\|_1$ (sparse).

Summary¶

• Regression: $A$ = samples × features, $\mathbf{x}$ = dense weights, $\mathbf{b}$ = targets. Focus on prediction.

• Sparse Representation: $A$ = signal × atoms, $\mathbf{x}$ = sparse activations, $\mathbf{b}$ = signal. Focus on compression/denoising.

• Shared Core: Both solve $A\mathbf{x} \approx \mathbf{b}$, but with different constraints on $\mathbf{x}$ and interpretations of $A$.

Connection Between Sparse Learning and Functional Analysis (Further Reading)¶

Sparse recovery techniques like OMP and LASSO are deeply rooted in functional analysis, which explains why terms like Hilbert spaces and Banach spaces frequently appear in machine learning research. These mathematical frameworks provide the theoretical foundation for understanding sparsity-inducing methods.

Hilbert Spaces and Orthogonal Projections¶

Hilbert spaces (complete inner product spaces) are essential for analyzing orthogonal projections and basis expansions. For example:

• OMP iteratively projects the residual onto the span of selected basis vectors, leveraging the orthogonality principle to minimize error at each step.

• The inner product structure of Hilbert spaces enables efficient computation of correlations between residuals and dictionary atoms (columns of $X$).

Banach Spaces and $\ell_1$-Regularization¶

LASSO operates in the context of Banach spaces (complete normed vector spaces) due to its reliance on the $\ell_1$-norm:

• The $\ell_1$-norm’s non-smoothness at the origin induces sparsity, a property studied in Banach space geometry.

• Unlike Hilbert spaces, Banach spaces generalize optimization techniques to non-Euclidean settings, crucial for sparse regularization.

Fourier, Wavelet, and Other Representations¶

Sparse learning often exploits transform domains (e.g., Fourier, wavelet, or learned dictionaries) where signals admit concise representations: - These domains provide structured bases where only a few coefficients are significant. - The connection to functional analysis arises because such bases typically form frames or Riesz bases in Hilbert spaces, ensuring stable sparse approximations.

From Fixed Bases to Learned Representations¶

1. Classical Sparse Coding:

   • Relies on predefined bases (e.g., Fourier, wavelets) where signals admit sparse representations.

   • These bases often form frames or Riesz bases in Hilbert spaces, ensuring stable approximations.

2. Convolutional Neural Networks (CNNs):

   • Learn adaptive bases/filters from data, implicitly constructing sparse-like representations through:

     • Local connectivity: Filters act as localized basis functions.

     • Activation sparsity: ReLU promotes de facto sparsity in feature maps.

   • While CNNs lack explicit $\ell_0$/$\ell_1$ constraints, their hierarchical structure approximates multi-scale sparse decompositions, akin to wavelet analysis but data-driven.

The Evolving Role of Sparse Methods¶

While deep learning has reduced reliance on handcrafted sparse models in some domains, sparse methods remain relevant because:

1. Interpretability:

   • OMP/LASSO yield explicit basis selections, whereas CNNs operate as black boxes.

   • Critical in fields like medicine or physics where model transparency is required.

2. Data-Efficiency:

   • Sparse methods often outperform DL in low-data regimes (e.g., medical imaging with small datasets).

3. Theoretical Guarantees:

   • Compressed sensing (OMP) and convex optimization (LASSO) provide recovery guarantees under precise conditions, unlike empirical DL results.

4. Hybrid Approaches:

   • Modern architectures (e.g., ISTA-Net, Learned Iterative Shrinkage) blend sparse priors with deep learning, showing that sparsity remains a useful inductive bias.

Sparse learning and deep learning are complementary:

• CNNs dominate when data is abundant and interpretability is secondary.

• OMP/LASSO persist in scenarios requiring rigor, efficiency, or transparency.

• Functional analysis bridges these paradigms, providing tools to analyze both fixed and learned representations.

MRI reconstruction (Further Reading)¶

MRI reconstruction can be formulated as an inverse problem where the goal is to recover the original image $\mathbf{m}$ from under-sampled k-space measurements $\mathbf{y}$. This is achieved using sparse optimization techniques.

A common formulation for MRI reconstruction is:

where:

• $\mathbf{m}$ represents the MRI image to be reconstructed.

• $\mathbf{F}_u$ is the under-sampled Fourier operator, mapping the image to k-space.

• $\mathbf{y}$ is the acquired k-space data (limited measurements).

• $\Psi$ is the sparsifying transform (e.g., wavelet or total variation).

• $\|\Psi \mathbf{m}\|_p$ represents the sparsity constraint, often chosen as $\ell_1$ (similar to LASSO) or $\ell_0$ (similar to OMP).

• $\lambda$ controls the balance between data fidelity and sparsity enforcement.

By leveraging sparsity, MRI reconstruction enables faster scans with fewer measurements while preserving essential anatomical details.

Source: Iterative reconstruction: how it works, how to apply it