F24 ECE 551

HW05, Due 11PM Thu. Oct. 10 1

Pr. 1.

Linear regression has a myriad of uses, including investigation of social justice issues.

Demo 05/sat-regress has data collected by the College Board–the organization that runs the SAT exam for high-school students, from Lesser, 2017. This data includes average SAT Math scores for 10 different family annual income brackets. This problem uses this data to explore the relationship between income and SAT scores.

Make a scatter plot of the midpoint of each income bracket versus the SAT Math score. The last income bracket says “100+” which does not have a midpoint, so just set it to be 120 (K$). Use the template code in the demo.

The scatter plot suggests there is a strong correlation between income and SAT score. You are going to fit four different models to the data (one at a time) and examine the fits and the coefficients:

For each of the four models, use a LLS fit to determine the β coefficient(s).

(a) Make a single plot showing the data and the four fits, for incomes between 0 and 130 (K$), i.e., income = 0:130 .

(We cannot predict anything for even higher incomes with this data.) Be sure to label your axes and use a legend to explain which points/lines are which.

(b) Make a table to report your β coefficients that looks like this:

(c) (Optional) In a statistics class, you would analyze the coefficients β1 and β2 to assess whether they are significantly different from 0 and establish evidence of a relationship. Even without formal statistics, you should be able to look at your plot and your table of coefficients and draw some conclusions about how equitable the SAT exam is.

Pr. 2.

Polynomial regression application

Let f(t) = 0.5 e 0.8t , t ∈ [0, 2].

(a) Suppose you are given 16 exact measurements of f(t), taken at the times t in the following 1D array:

T = range(0, 2, 16)

Use Julia to generate 16 exact measurements: bi = f(ti), i = 1, . . . , 16, for ti ∈ T .

Now determine the coefficients of the least-squares polynomial approximation of the data b for

(i) a polynomial of degree 15: p15(t);

(ii) a polynomial of degree 2: p2(t).

Compare the quality of the two approximations graphically. Use scatter to first show bi vs ti for i = 1, . . . , 16, then use plot! to add plots of p15(t) and p2(t) and f(t) to see how well they approximate the function on the interval [0, 2]. Pick a very fine grid for the interval, e.g., t = (0:1000)/500 . Make the y-axis range equal [−1, 4] by using the ylim=(-1,4) option. As always, include axis labels and a clear legend. Submit your plot and also summarize the results qualitatively in one or two sentences.

For this problem, write your own Julia code rather than using fit in the Polynomials package. (You may check your solutions using that function.)

(b) Now suppose the measurements are affected by some noise. Generate the measurements using yi = f(ti) + ei , i = 1, . . . , 16, , as follows. You must use the seed to get the correct values!

using Random: seed!

seed!(3); e = randn(length(T))

Determine the coefficients of the least-squares polynomial approximation of the (noisy) measurements y for

(i) a polynomial of degree 15: p15(t);

(ii) a polynomial of degree 2: p2(t).

Compare the two approximations as in part (a). Again make the y-axis range equal [−1, 4] by using the ylim=(-1,4) option. Submit your plot and also summarize the results qualitatively in a couple of sentences, including comparing the behavior. in (a) and (b).

(c) Let xˆn(b) and xˆn(y) denote the LLS polynomial coeffients from noiseless b and noisy y, respectively, for a polynomial of degree n. Report the values of the residual norms ∥Axˆn(b) − b∥2 and ∥Axˆn(y) − y∥2 for the polynomial fits of degree 2 and degree 15. These residual norms describe how closely the fitted curve fits the data.

Also report the fitting errors ∥Axˆn(y) − b∥2 for n = 2, 15. Arrange your results in a table as follows:

Hint: the bottom left entry is between 1.6 and 1.9.

(d) Explain why the residual norm for degree 2 is smaller or larger than that of degree 15.

Pr. 3.

(a) One use of the “vec trick” vec(AXB⊤) = (B ⊗ A)vec(X) is computing a 2D discrete Fourier transform. (DFT) of a 2D signal. The (1D) DFT of a vector x ∈ C N is the vector fx ∈ C N with entries

(Here we follow linear algebra (and Julia) where the first array index is 1, whereas DSP books usually use 0, . . . , N − 1.) We can represent the above computation in matrix-vector form. as fx = FN x, where FN is the N × N DFT matrix with entries

We can generate FN in Julia as follows.

# N × N DFT matrix

using FFTW: fft

using LinearAlgebra: I

F = fft(I(N), 1)

One can verify that FN ′ FN = NIN , so the (1D) inverse DFT of fx can be computed as x = (1/N)FN ′ fx.

We compute the 2D DFT, call it SX, of the M × N matrix X by computing the 1D DFT of each column of X followed by the 1D DFT of each row of the result (or vice versa).

Explain why the following expression computes the 2D DFT of X:

(b) Write vec(SX) as the product of a matrix and vec(X).

(c) Show that the following expression computes the 2D inverse DFT of SX:

where Y denotes (elementwise) complex conjugate of matrix Y .

Hint: Use the fact that FN ′ FN = FN FN ′ = NIN .

(d) Write vec(X) as the product of a matrix and vec(SX).

Pr. 4.

In the ordinary least-squares problem we found the (minimum norm) xˆ that minimized ∥r∥2 where r = b − Ax is the residual. We saw that the optimal x can be expressed entirely in terms of b and an SVD of A. Let A be an M × N matrix so that x ∈ F N and b ∈ FM and r ∈ FM.

In applications with heteroscedastic measurement errors, we prefer to minimize ∥W r∥2, i.e., a weighted squared error, where W is a diagonal matrix having diagonal entries w1, . . . wM ≥ 0. Determine the optimal x for this weighted least-squares problem. (Again, you should express the answer in terms of b, W, and an SVD of an appropriate matrix.) Your answer must not have any pseudoinverse in it! (It may have an inverse as long as you are sure that the matrix is invertible and trivial to invert.) You may assume that W*A == zeros(size(W*A)) is false in Julia.

Pr. 5.

(a) Find (by hand) all solutions of the linear least-squares problem arg minx∈R2 ∥Ax − b∥2 when  and 

(b) Describe briefly how you would modify your solution for this variation of the problem: arg minx∈C2∥Ax − b∥2.

Pr. 6.

Recall that the least-squares problem

has the solution xˆ = A+b where A+ denotes the pseudoinverse. When A is large, it can be expensive to compute xˆ using the pseudoinverse of A. In such settings, the gradient descent (GD) iteration given by

xk+1 = xk − αA′ (Axk − b),                                    (1)

will converge to a minimizer of ∥Ax − b∥2 as iteration k → ∞ when 0 < α < 2/σ1 2 (A). Note that the iteration has a fixed point, i.e., xk+1 = xk if

A′ (Axk − b) = 0,

which are exactly the normal equations. So any fixed point minimizes the least-squares cost function.

(a) Write a function called lsgd that implements the above least-squares gradient descent algorithm.

In Julia, your file should be named lsgd.jl and should contain the following function:

"""

x = lsgd(A, b ; alpha=0, x0=zeros(size(A,2)), nIters::Int=200)

Perform. gradient descent to solve the least-squares problem:

``\\argmin_x 0.5 ‖ A x - b ‖2``

# In:

- `A` `m × n` matrix

- `b` vector of length `m`

# Option:

- `alpha` step size to use, and must satisfy ``0 < α < 2 / σ1(A)^2``

to guarantee convergence,

where ``σ1(A)`` is the first (largest) singular value.

Ch.6 explains a default value for `alpha`

- `x0` is the initial starting vector (of length `n`) to use.

Its default value is all zeros for simplicity.

- `nIters` is the number of iterations to perform. (default 200)

Out:

- `x` vector of length `n` containing the approximate LS solution

"""

function lsgd(A, b ; alpha::Real=0, x0=zeros(size(A,2)), nIters::Int=200)

Email your solution as an attachment to [email protected].

The function specification above uses a powerful feature of Julia where functions can have optional arguments with specified default values. If you simply call lsgd(A,b) then alpha , x0 and nIters will all have their default values. But if you call, say, lsgd(A, b, nIters=5, alpha=7) then it will use the specified values for nIters and alpha and the default for x0 . Note that these named optional arguments can appear in any order. This approach is very convenient for functions with multiple arguments.

(b) Use your function (after it passes) to generate a plot of log10(∥xk − xˆ∥) versus k = 0, 1. . . . , 200 using α = 1/σ1 2 (A) for A and b generated as follows.

using Random: seed!

m, n = 100, 50; sigma = 0.1

seed!(0) # seed random number generator

A = randn(m, n); xtrue = rand(n); noise = randn(m)

b = A * xtrue + sigma * noise # b and xhat change with σ

Repeat for σ = 0.5, 1, 2 and submit one plot with all four curves on it. Does ∥xk − xˆ∥ decrease monotonically with k in the plots? Note that σ = sigma here is a noise standard deviation unrelated to singular values.

Pr. 7.

(a) For δ > 0, determine the solution of the regularized LS problem

Express the answer in terms of the SVD of an appropriate matrix.

Use the fact that B+ = (B′B) −1B′ when B′B is invertible to simplify the expression.

Hint: Rewrite the cost function so it looks like a usual least-squares problem.

(b) What does xˆ(δ) tend to as δ → ∞? Does this answer make sense?

(c) Write (by hand, not code) an iteration based on the gradient descent (GD) method such that the iterates converge to the minimizer xˆ(δ).

(d) Determine a condition on the step size µ that guarantees convergence of your GD method. Express the condition in terms of the original problem quantities: A, b, and δ.

Non-graded problem(s) below

(Solutions will be provided for self check; do not submit to gradescope.)

Pr. 8.

Let S and T denote subspaces of a vector space V. Prove or disprove the following.

(a) If S is a subset of T or vice versa, then the union of subspaces S ∪ T is a subspace of V.

(b) If the union of subspaces S ∪ T is a subspace of V, then S is a subset of T or vice versa.