Department of Electrical Engineering

Indian Institute of Technology Kharagpur

Algorithms and AI-ML Laboratory (EE22202)

Assignment 2: Regression

Krish Agrawal
Roll No: 24IE10027

February 12, 2026

Contents

2
Question 1: Feature Map Definition
Problem: Define a suitable feature map ϕ(x) which maps x ∈ R4 to entries of a polynomial
of degree 2 in x. What is the dimension of ϕ(x)?

Solution
Let the input vector be x = [x1 , x2 , x3 , x4 ]⊤ . A general polynomial of degree 2 includes the
constant term (degree 0), linear terms (degree 1), and quadratic terms (degree 2).
The feature map ϕ(x) collects all these terms:

1. Degree 0 (Bias): The constant term 1. (1 term)

2. Degree 1 (Linear): The variables x1 , x2 , x3 , x4 . (4 terms)

3. Degree 2 (Quadratic):

• Squared terms: x21 , x22 , x23 , x24 . (4 terms)

• Cross-products: x1 x2 , x1 x3 , x1 x4 , x2 x3 , x2 x4 , x3 x4 . (6 terms)

Total dimension calculation:

Dimension = 1 + 4 + 4 + 6 = 15

Alternatively, using the binomial coefficient for n = 4 variables and degree d = 2:

     
n+d 4+2 6
= = = 15
d 2 2

Answer: The dimension of ϕ(x) is 15.

Question 2: Linear Regression Formulation

Problem: Formulate a linear regression problem to determine the coefficients of this polyno-
mial. Clearly state the decision variable w, its dimension, and the cost function.

Solution
We formulate the problem as a linear regression on the transformed features.

1. Decision Variable
The decision variable is the weight vector w, which contains a coefficient for every term in the
feature map ϕ(x). Since ϕ(x) has dimension 15:

w ∈ R15

1
2. Cost Function
We assume a dataset of N = 100 samples given by pairs {x̂i , ŷi }. We define the design matrix
Φ ∈ R100×15 where the i-th row is the feature map ϕ(x̂i ).
The objective is to minimize the sum of squared errors (Least Squares):
100
X
min15 L(w) = (ŷi − ϕ(x̂i )w)2
w∈R
i=1

In vector notation:
min15 ∥ŷ − Φw∥22
w∈R

Question 3
Problem: Compute the gradient of this cost function with respect to the decision variable w.

Solution
Let the cost function be denoted by l(w). From Question 2, we have the vector form:

l(w) = ∥ŷ − Φw∥22

Using the vector identity ∥v∥22 = v ⊤ v, we expand the objective function:

l(w) = (ŷ − Φw)⊤ (ŷ − Φw)

= (ŷ ⊤ − w⊤ Φ⊤ )(ŷ − Φw)
= ŷ ⊤ ŷ − ŷ ⊤ Φw − w⊤ Φ⊤ ŷ + w⊤ Φ⊤ Φw

Note that ŷ ⊤ Φw is a scalar quantity. Since the transpose of a scalar is the scalar itself, we have
(ŷ ⊤ Φw)⊤ = w⊤ Φ⊤ ŷ. Thus, we can combine the middle terms:

l(w) = ŷ ⊤ ŷ − 2w⊤ Φ⊤ ŷ + w⊤ Φ⊤ Φw

Now, we compute the gradient ∇w l(w) by differentiating each term with respect to w:
1. The derivative of the constant term ŷ ⊤ ŷ is 0.

2. The derivative of the linear term −2w⊤ (Φ⊤ ŷ) is −2Φ⊤ ŷ (using ∇x (x⊤ a) = a).

3. The derivative of the quadratic term w⊤ (Φ⊤ Φ)w is 2Φ⊤ Φw (using ∇x (x⊤ Ax) = 2Ax for
symmetric matrix A = Φ⊤ Φ).
Combining these results:
∇w l(w) = −2Φ⊤ ŷ + 2Φ⊤ Φw
Factoring out 2Φ⊤ , we obtain the final gradient expression:

∇l(w) = 2Φ⊤ (Φw − ŷ)

Question 4: Optimal Weights and Polynomial

Problem: Find the optimal weights w∗ using the Normal Equation and determine the poly-
nomial that maps input to output.

2
 Methodology
Since the objective function l(w) = ∥ŷ − Φw∥22 is a convex quadratic function, the optimal
weights w∗ can be found analytically by setting the gradient to zero. This yields the Normal
Equation:
w∗ = (Φ⊤ Φ)−1 Φ⊤ ŷ
This formula provides the exact solution for the unconstrained least squares problem without
requiring an iterative QP solver.

Python Implementation
The following code first generates the synthetic dataset (as per the problem description) and
then solves for w∗ using the Normal Equation.
Listing 1: Data Generation and Solution via Normal Equation
1 import numpy as np
2

3 # --- 1. Data Generation ---

4 np . random . seed (42) # Fixed seed for reproducibility
5

6 N = 100 # Number of samples

7 d = 4 # Input dimension
8

9 # Generate random input X ( N x 4)

10 X_data = np . random . randn (N , d )
11

12 # True hidden weights ( for generating target y )

13 w_true = np . random . randn (15) * 0.8
14

15 def get_feature_map ( x ) :
16 Maps 4 D input to 15 D polynomial feature vector ( Degree 2)
17 x0 , x1 , x2 , x3 = x
18 return np . array ([
19 1, # Bias
20 x0 , x1 , x2 , x3 , # Linear
21 x0 **2 , x1 **2 , x2 **2 , x3 **2 , # Squared
22 x0 * x1 , x0 * x2 , x0 * x3 , # Cross - product
23 x1 * x2 , x1 * x3 ,
24 x2 * x3
25 ])
26

27 # Create Design Matrix Phi ( N x 15)

28 Phi = np . vstack ([ get_feature_map ( x ) for x in X_data ])
29

30 # Generate Target y = Phi * w_true + noise

31 noise = 0.1 * np . random . randn ( N )
32 y = Phi @ w_true + noise
33

34 # --- 2. Solve using Normal Equation ---

35 # Formula : w * = ( Phi ^ T * Phi ) ^( -1) * ( Phi ^ T * y )
36

37 XtX = Phi . T @ Phi

38 Xty = Phi . T @ y

3
39

40 # np . linalg . solve is numerically more stable than inv ()

41 w_star = np . linalg . solve ( XtX , Xty )
42

43 print ( Optimal Weights w *: , np . round ( w_star , 4) )

Resulting Polynomial
The determined polynomial f (x) is given by the dot product of the optimal weight vector and
the feature map:
X15
∗
f (x) = ϕ(x)w = wj∗ ϕj (x)
j=1

Substituting the learned weights, the function approximates the underlying relationship be-
tween the input variables x1 , . . . , x4 and the output y.

Question 5: Error Vector and Histogram

Problem: Compute the error vector ŷ − Φw∗ and plot its histogram.

Methodology
The error vector (or residual vector) e ∈ RN is defined as the difference between the observed
target values and the values predicted by our model:

e = ŷ − Φw∗

Since the synthetic data was generated as y = Φwtrue + noise, where the noise follows a normal
distribution N (0, σ 2 ), we expect the histogram of the residuals to approximate a Gaussian
distribution centered around 0.

Python Implementation
The following code computes the residuals and plots the histogram using Matplotlib.
Listing 2: Computing Residuals and Plotting Histogram
1 import matplotlib . pyplot as plt
2

3 # 1. Compute the Error Vector ( Residuals )

4 # residuals = y_actual - y_predicted
5 residuals = y - Phi @ w_star
6

7 # 2. Plot the Histogram

8 plt . figure ( figsize =(8 , 6) )
9 plt . hist ( residuals , bins =20 , color = ' skyblue ' , edgecolor = ' black ' , alpha
,→ =0.7)
10

11 # Labels and Title

12 plt . xlabel ( ' Residual Error ( $y - \ hat { y } $ ) ')
13 plt . ylabel ( ' Frequency ')
14 plt . title ( ' Histogram of Regression Residuals ')

4
15 plt . grid ( True , linestyle = ' -- ' , alpha =0.5)
16

17 # Save the plot

18 plt . savefig ( ' histogram . png ')
19 plt . show ()

Resulting Plot
The histogram below visualizes the distribution of the errors.

Figure 1: Histogram of the residual errors. The distribution is centered near zero, consistent
with the Gaussian noise added during data generation.

Question 6: Optimization Algorithms (GD, AGD, SGD)

Problem: Solve the regression problem using Gradient Descent (GD), Accelerated Gradient
Descent (AGD), and Stochastic Gradient Descent (SGD) for 1000 iterations. Plot the log of
the cost function and the log of the error norm ∥wt − w∗ ∥ versus iterations.

5
 Methodology
We initialized all algorithms at w0 = 0. The optimal solution w∗ calculated in Question 4 is
used as the ground truth for error tracking.
Step Size Selection:

• GD and AGD: The step size η was set based on the Lipschitz constant L of the gradient
(maximum eigenvalue of the Hessian H = 2Φ⊤ Φ). We used η = L1 .

• SGD: Since SGD approximates the gradient using a single sample, it requires a dimin-
η0
ishing step size to converge. We used a schedule ηt = 1+αt with η0 = 0.01.

Python Implementation
The following code implements the three algorithms and generates the convergence plots.
Listing 3: Implementation of GD
1 import numpy as np
2 import matplotlib . pyplot as plt
3

4 # ... ( Data and w_star from previous steps ) ...

5

6 # 1. Setup Parameters
7 # Hessian H = 2 * Phi . T @ Phi
8 eigenvalues = np . linalg . eigvalsh (2 * Phi . T @ Phi )
9 L = np . max ( eigenvalues ) # Lipschitz constant
10

11 eta_gd = 1.0 / L
12 eta_agd = 1.0 / L
13 eta_sgd = 0.01
14 T = 1000
15 w_init = np . zeros (15)
16

17 # 2. Gradient Descent
18 w = w_init . copy ()
19 gd_loss , gd_dist = [] , []
20 for t in range ( T ) :
21 grad = 2 * Phi . T @ ( Phi @ w - y )
22 w = w - eta_gd * grad
23 gd_loss . append ( np . sum (( y - Phi @ w ) **2) )
24 gd_dist . append ( np . linalg . norm ( w - w_star ) )
25

26 # 3. Accelerated Gradient Descent ( Nesterov )

27 w = w_init . copy ()
28 y_nest = w_init . copy ()
29 agd_loss , agd_dist = [] , []
30 for t in range ( T ) :
31 # Gradient step
32 grad = 2 * Phi . T @ ( Phi @ y_nest - y )
33 w_next = y_nest - eta_agd * grad
34

35 # Momentum step
36 beta = t / ( t + 3)
37 y_nest = w_next + beta * ( w_next - w )

6
38

39 w = w_next
40 agd_loss . append ( np . sum (( y - Phi @ w ) **2) )
41 agd_dist . append ( np . linalg . norm ( w - w_star ) )
42

43 # 4. Stochastic Gradient Descent

44 w = w_init . copy ()
45 sgd_loss , sgd_dist = [] , []
46 for t in range ( T ) :
47 i = np . random . randint (0 , N )
48 grad_i = 2 * ( np . dot ( Phi [ i ] , w ) - y [ i ]) * Phi [ i ]
49

50 # Diminishing learning rate

51 lr = eta_sgd / (1 + 0.01 * t )
52 w = w - lr * grad_i
53

54 sgd_loss . append ( np . sum (( y - Phi @ w ) **2) )

55 sgd_dist . append ( np . linalg . norm ( w - w_star ) )
56

57 # Plotting code omitted for brevity ( see Resulting Plots )

Resulting Plots
The figures below compare the convergence behavior.

Figure 2: Left: Logarithm of the cost function vs. iterations. Right: Logarithm of the weight
error ∥wt − w∗ ∥ vs. iterations.

Observations
• AGD (Accelerated Gradient Descent): Converges the fastest, reaching a low er-
ror significantly earlier than standard GD. The momentum term effectively accelerates
progress along shallow curvature directions.

• GD (Gradient Descent): Shows steady, linear convergence (straight line on the log
plot) but is slower than AGD.

7
 • SGD (Stochastic Gradient Descent): Descends rapidly in the very early iterations
but exhibits ”noisy” behavior. The curve fluctuates because the gradient is estimated
from single samples, making it much rougher than the full-batch methods.

Question 7: Convergence vs. Gradient Evaluations

Problem: Plot the (log of) cost function and weight error against the number of gradient
evaluations.

Methodology
To compare the computational efficiency, we transform the x-axis from ”iterations” to ”gradient
evaluations”:

• GD and AGD: Compute the gradient on the full dataset (N = 100) at each step.
Therefore, at iteration t, the number of evaluations is t × 100.

• SGD: Computes the gradient on a single sample at each step. Therefore, at iteration t,
the number of evaluations is t × 1.

Python Implementation
We adapt the plotting logic to use these transformed x-axis values.
Listing 4: Plotting against Gradient Evaluations
1 # Transformation Logic
2 # GD / AGD : x - axis = iteration * 100
3 # SGD : x - axis = iteration * 1
4

5 plt . figure ( figsize =(14 , 6) )

6

7 # Left : Cost vs Evaluations

8 plt . subplot (1 , 2 , 1)
9 plt . semilogy ( hist_gd [ ' evals '] , hist_gd [ ' cost '] , label = ' GD ')
10 plt . semilogy ( hist_agd [ ' evals '] , hist_agd [ ' cost '] , label = ' AGD ')
11 plt . semilogy ( hist_sgd [ ' evals '] , hist_sgd [ ' cost '] , label = ' SGD ')
12 plt . xlabel ( ' Gradient Evaluations ')
13 plt . ylabel ( ' Cost ( Log Scale ) ')
14 plt . legend ()
15

16 # Right : Error vs Evaluations

17 plt . subplot (1 , 2 , 2)
18 plt . semilogy ( hist_gd [ ' evals '] , hist_gd [ ' dist '] , label = ' GD ')
19 plt . semilogy ( hist_agd [ ' evals '] , hist_agd [ ' dist '] , label = ' AGD ')
20 plt . semilogy ( hist_sgd [ ' evals '] , hist_sgd [ ' dist '] , label = ' SGD ')
21 plt . xlabel ( ' Gradient Evaluations ')
22 plt . ylabel ( ' Weight Error ( Log Scale ) ')
23 plt . legend ()
24

25 plt . show ()

8
 Resulting Plots

Figure 3: Convergence efficiency. Note that for the same number of iterations (1000), GD
and AGD perform 100,000 evaluations (extending far right), while SGD performs only 1,000
(compressed on the left). This illustrates the high computational cost per step of full-batch
methods.

Question 8: Under-determined System (N = 10)

Problem: Compute optimal weights w∗ using only the first 10 data points. Plot the error
histogram and discuss the uniqueness of the solution.

Methodology
In this case, we have N = 10 samples and D = 15 features (weights). Since N < D, the system
of equations Φw = y is under-determined. The matrix Φ⊤ Φ is of size 15 × 15 but has a
maximum rank of 10. Therefore, it is singular (non-invertible) and has a non-trivial null space.
Uniqueness: The optimal solution w∗ is not unique. There is an affine subspace of
solutions of dimension at least D − N = 5. While a QP solver or pseudo-inverse will return a
single specific solution (typically the one with the minimum L2 norm, ∥w∥2 ), infinitely many
other solutions exist that achieve zero training error.

Python Implementation
We subset the data and solve using a least-squares solver (which handles singular matrices).
Listing 5: Solving the Under-determined System
1 import numpy as np
2 import matplotlib . pyplot as plt
3

4 # --- 1. Subset Data ( First 10 points ) ---

5 Phi_sub = Phi [:10] # Shape (10 , 15)
6 y_sub = y [:10] # Shape (10 ,)
7

8 # --- 2. Solve for w * ---

9
 9 # Since Phi is not full rank , we use lstsq
10 # ( finds min - L2 - norm solution among infinite solutions )
11 w_sub , residuals , rank , s = np . linalg . lstsq ( Phi_sub , y_sub , rcond = None
,→ )
12

13 print ( f Rank of Design Matrix : { rank } )

14 print ( f System is Under - determined : { rank } < 15 )
15

16 # --- 3. Plot Error Histogram ( Training Error ) ---

17 # For N < D , we expect training error to be practically zero (
,→ overfitting )
18 train_error = y_sub - Phi_sub @ w_sub
19

20 plt . figure ( figsize =(8 , 6) )

21 plt . hist ( train_error , bins =5 , color = ' orange ' , edgecolor = ' black ')
22 plt . title ( ' Error Histogram ( N =10 Subset ) ')
23 plt . xlabel ( ' Residual Error ')
24 plt . ylabel ( ' Frequency ')
25 plt . grid ( True , alpha =0.3)
26 plt . show ()

Question 9: L1 Regularization (Lasso)

∗
Problem: Modify the cost function by adding an L1 regularization term and solve for wreg
using CVXPY for different values of λ.

Formulation
The Lasso (Least Absolute Shrinkage and Selection Operator) optimization problem is:

min ∥y − Φw∥22 + λ∥w∥1

w
P
Unlike the L2 norm, the L1 norm (∥w∥1 = |wi |) is not differentiable at zero, so we use a
convex solver like CVXPY.

Python Implementation

Listing 6: Lasso Regression with CVXPY

1 import cvxpy as cp
2

3 lambdas = [0.1 , 1.0 , 10.0]

4 results = {}
5

6 print ( f { ' Lambda ': <10} | { ' Sparsity ( Zeros ) ': <15} | { ' L2 Error ': <10} )
7 print ( - * 45)
8

9 for lam in lambdas :

10 # Define Variables
11 w_var = cp . Variable (15)
12

13 # Define Objective : MSE + Lambda * L1_norm

10
14 loss = cp . sum_squares ( y_sub - Phi_sub @ w_var )
15 reg = lam * cp . norm ( w_var , 1)
16 objective = cp . Minimize ( loss + reg )
17

18 # Solve
19 prob = cp . Problem ( objective )
20 prob . solve ()
21

22 w_lasso = w_var . value

23

24 # Count zeros ( sparsity )

25 # We use a small threshold because floating point math is rarely
,→ exactly 0
26 n_zeros = np . sum ( np . abs ( w_lasso ) < 1e -4)
27 results [ lam ] = w_lasso
28

29 print ( f { lam : <10} | { n_zeros : <15} | { prob . value :.4 f } )

30

31 # Example : Inspecting weights for lambda =1.0

32 print ( \ nLasso Weights ( lambda =1.0) : )
33 print ( np . round ( results [1.0] , 3) )

Question 10: Comparison of Solutions

Problem: Comment on the difference between the solutions w∗ (from Q8) and wreg
∗
(from Q9).

Discussion
Comparing the unregularized solution (w∗ ) and the L1-regularized solution (wreg
∗
):

1. Sparsity:
• w∗ (Minimum Norm Solution): This solution is typically dense. It tends to dis-
tribute small non-zero values across all 15 weights to minimize the L2 norm while
satisfying the training data perfectly.
• wreg
∗
(Lasso Solution): This solution is sparse. The geometry of the L1 penalty (a
diamond shape) encourages the optimization solution to land on the axes, effectively
setting irrelevant feature weights to exactly zero.
2. Interpretability:
• wreg
∗
performs feature selection, identifying the polynomial terms that are most
critical for mapping x to y. This makes the Lasso model easier to interpret.
• w∗ uses all features, making it harder to distinguish signal from noise.
3. Overfitting:
• w∗ perfectly memorizes the 10 training points (training error ≈ 0) but likely gener-
alizes poorly to new data.
• wreg
∗
sacrifices some training accuracy (slightly higher training error) to keep the
weights simple, often leading to better generalization on unseen data.

11