CAAM 454/554 Homework Assignments

The problems are numbered as in the textbooks. TB refers to Trefethen and Bau, and NW to Nocedal and Wright (2nd Ed.).

Problems starting with (554) will be required only from CAAM 554 students (similarly for those with (454)).

HW 10: Due 4/11, Thur. (total points: 454/554 = 50/65) (Solution)

(40 points) Assignment: Primal-Dual Interior-point method for LP
Submit your output (truncated the iteration history if too long) along with your code and a short write-up (typed) explaining what you observe.
(5 points each) Problems 14.8, 14.9.
(554, 5 points each) Problems NW 13.5, 14.11, 14.17.

HW 9: Due 4/2, Tue. (total points: 454/554 = 30/30) (Solution)

Consider again the low-rank approximation nonlinear least squares problem.
This time the real Gauss-Newton method without line search:
```
      Y = X/(X'*X); 
      P = (I - 0.5*X*Y')*(A*Y - X);
      X = X + P;
```
--- (10 pts) Write a function to implement the above scheme
```
       [X,iter] = myGN4real(A,k,tol,maxit)
```
that solves the NLS problem. Optimize your code as much as you can. Stop the iterations when relchg = norm(P,'fro)/norm(X,'fro') is less than tol. Mimic the on-screen output format of the instructor's code.
--- (10 pts) Write a function to evaluate the objective function. Try to use computed quantities (by your main function) most efficiently (about 2k^2 scalar multiplications would be necessary).
--- Download and run the test script testGN4real.m with the instructor's code yzGN4real.p (output should look like output.pdf for k = (1:10)*50). Also run with quiet = false to see on-screen printout.
--- Edit the script to run your code. First try small k (say, k = (1:5)*20). Submit your code along with (a) the output of k = (1:10)*50 for quiet = true, and (b) the output of k = 600 for quiet = false (or a reduced k value if necessary).
--- (10 pts) Submit a short write-up (typed) explaining what you observe. Use the Matlab profiler to find out and then report the computational cost profile of your code in the runs (a) and (b) above.

HW 8: Due 3/21, Thur. (total points: 454/554 = 55/75) (Solution)

(20+5 pts) Write a matlab function for BFGS method with the following interface:
```
function [x,iter,hist] = myBFGS(func,x0,tol,maxit,varargin)
%
% Quasi-Newton method using inverse BFGS update 
% (if you do backtracking, make sure you check y'*s > 0)
%
% Usage: 
%   [x,iter,hist] = myBFGS(func,x0,tol,maxit,varargin);
%
% INPUT:
%       func - matlab function for f and g evaluations
%         x0 - starting point (column vector)
%        tol - stopping tolerance for gradient norm
%      maxit - maximum number of iterations
%
% OUTPUT:
%          x - computed approximate solution
%       iter - number of iterations taken
%      hist(1,:) - vector of function values at all iterations
%      hist(2,:) - vector of gradient norms  at all iterations
```
--- Use a multiple of identity matrix as the initial inverse Hessian approximation.
--- Stop when the norm of gradient is less than the tolerance times (1+abs(f)).
--- Make your output format similar with that of the instructor's code (for grading purposes).
--- Download the instructor's code yzBFGS.p and yzCGD.p (then you can run the test scripts).
--- Run the script testBFGS.m, with your code myBFGS on the function rosenbrock_g.m (generalized Rosenbrock).
--- Run the script testDenoise2.m with your codes myBFGS, myCDG on the function TVL2N.m (from a previous assignment).
--- (5 pts) Submit the output (truncated printout and graphs) along with your code and a short write-up (typed) explaining what you observe.
(15 pts) Solve the low-rank approximation problem: for given n by n, symmetric positive semi-definite A and k < n, find
```
      X = argmin {0.25||XX^T - A||_F²/||A||_F²: X is n by k}
```
--- Write a function [f,g] = mylork(x,n,k,A,nrmA2) where the small x is the vectorized big X and nrmA2 = norm(A,'fro')^2.
--- Run the script testBFGS2.m, with your code myBFGS on the function.
--- Compare with the instructor's code yzlork.p and yzBFGS.p (can be run without your codes).
--- Submit the output along with your code and a short write-up (typed) explaining what you observe.
(15 pts) The above low-rank approximation problem is a nonlinear least squares problem.
When X has full rank, a Gauss-Newton like method with line search takes the form: choose a step size a, compute
```
      G = X*(X'*X) - A*X
      X = X - a * G / (X'*X)
```
--- Write a function to implement the above scheme
```
       [X,iter] = myGNes(A,k,tol,maxit)
```
that solves the NLS problem (you can assume norm(A,'fro')=1).
--- Download the image files Panda1.mat, Panda2.mat, and the instructor's code yzGNes.p.
--- Run the script testPanda.m with your code myGNes.m. Edit the script to run on the 2nd image.
--- Submit the output along with your code and a short write-up (typed) explaining what you observe.
(554, 20 pts) Problems NW 6.1(a), 6.3, 6.4, 6.9 (prove Psi(B)>= n instead).

HW 7: Due 3/14, Thur. (total points: 454/554 = 60/100) (Solution)

(40 pts) Problems NW 10.1 - 10.4.
(554, 20 pts) Problems NW 10.5 - 10.6.
(554 10 pts) Study Theorems 3.5, 11.2 and 11.3, understand what they say, their proofs and connections between them. Then stete and prove a theorem, following the model of Theorem 11.3, for the damped Newton method that solves the nonlinear equations r(x)=0. Specifically, the update using p that satisfies (11.17) is replaced by
```
x_k+1 = x_k + (1 - η_k)p_k, for some η_k in [0,η], 0 < η < 1,
```
where p_k is the Newton direction at x_k.
(554 10 pts) Consider the following algorithm for solving nonlinear systems. At the current iterate x_k, do a Newton iteration to generate an intermediate point. Then do another Newton-like iteration from this intermediate point but using the same Jacobian (evaluated at x_k) to get the next iterate. (a) Write the iteration into a fixed-point iteration form x_k+1 = S(x_k). (b) What is the rate of convergence for this algorithm (under suitable conditions)? Sketch a formal argument for your claim. (This part will not be graded. Just do your best.)
(20 points) Implement the Gauss-Newton method with a line search to solve a set of nonlinear least squares problems. The interface of the code should be:
```
 
    [x,iter] = MyGN(func,x,tol,maxiter,prt,varargin);
```
where func is a function name that returns [r,J], prt is a switch for printing iteration info or not, and varargin is a variable list of input arguments to be used by function func. The stopping criterion should be norm(g) < tol, where g = J'*r is the gradient vector.
Download and unzip the set of codes in nls.zip. Download the instructor's codes ZyGN.p and ZyLM.p. Study and run the test code testNLS.m. Without input, it runs the instructor's codes (CodeRange = 3:6) by default. We included 33 test problems, as described in the paper MGH.pdf. The test problems have already been coded. (A program FdNewton.m is included that implements Newton's method for solving nonlinear systems using a finite-difference scheme to approximate Jacobian. You may start your programming by studying and modifying this program.)
Submit your code and printout from the test program testNLS.m for the following two runs:
```
        testNLS(1, 1.e-4,[1 3]);   (easy run) 
        testNLS(10,1.e-6,[1 3]);   (hard run)
```

HW 6: 3/7 Thur. (total points: 454/554 = 70/100) (Solution)

All functions are assumed to be from Rⁿ to R unless otherwise specified.

(1) (10 pts) Show that the intersection of a collection of convex sets is convex.
(2) (10 pts) Show that if f(x) is a convex function, then the set {x: f(x) &le 0} is a convex set.
(3) (10 pts) Prove that the set of real symmetric positive semidefinite matrices is convex.
(4) (554 10 pts) Prove that the point-wise maximum of a (finite or infinite) set of convex functions is convex.
(5) (554 10 pts) Prove that a C² function is convex if and only if its Hessian is positive semidefinite at all x.
(6) (554 10 pts) Prove that the maximum eigenvalue is a convex function over the set of real symmetric positive semidefinite matrices.
(20 pts) Write a gradient (steepest) descent Matlab function:
```
	function [x,iter,nf] = myCGD(func,method,x0,tol,maxit,varargin)
```
(where method is a dummy parameter for now) to minimize a function defined by the function handle func starting from initial guess x0. Use back-tracking line search and choose your own parameters. The outputs include an approximate solution, the iteration number taken, and the number of function evaluations taken. You should mimic the printout format as in the instructor's code.
You interface func by
```
	[f,g] = feval(func,x,varargin{:}); % function and gradient
	f = feval(func,x,varargin{:});     % function value only
```
The stopping criterion is either
```
	crit := norm(g)/(1 + abs(f)) < tol;
```
or the maximum iteration number maxit is reached.
Download the test scripts testDenoise.m and testDeblur.m, the functions TVL2N.m and TVL2B.m, and the instructor's code yzCGD.p. You can run the test scripts without your codes in place.
Run the test scripts with your code (you will need image processing toolbox). Submit your code and printout including graphs, and a short report (typed) explaining what you observe in these experiments.
(10+10 pts) (a) Write a Matlab function to minimize C² convex functions using Newton's method (without line search):
```
	function [x,iter,nf] = myNewton(func,x0,tol,maxit,varargin).
```
Use the same stopping rule as above. The test function is a quartic f(x) = (x'*A*x)^2/4 + a*(sum(x) - n)^2/2.
You need to write a Matlab function [f,g,H] = Aquartic(x,A,a) to calculate the function, gradient and Hessian values of this function.
(b) Modify your myCGD.m code to also allow the use of nonlinear CG methods, where method = 1,2,3,4 will correspond to GD, CG-PR, CG-PR+ and CG-FR, respectively. If a search direction is non-descent, then use the steepest descent direction instead.
Download the instructior's code yzNewton.p and yzAquartic.p, and the test script testAquartic.m, which can be run without your code. When your codes are ready, run the test script for n = 1000.
Submit your codes and the printout including graphs, and a short write-up (typed) explaining what you observe in these experiments. Submit a short write-up (typed) explaining your observations in these experiments, such as convergence rates, with printout and other necessary supporting evidence.

HW 5: 2/21 Thur. (total points: 454/554 = 40/40) (Solution)

(10 pts) Write a Matlab function function X = mySolve4X(A,B,tol,maxit) to solve linear systems in a matrix form AX + XA = B for unknown matrix X, where all involved are square matrices. Assume that A is symmetric positive definite, and use the Matlab pcg function to solve the system by supplying a function (handle) to perform matrix-vector multiplications corresponding to the linear operation AX + XA. Turn off the on-screen printout of pcg function by specifying 2 outputs (see help pcg). You may use Matlab function reshape to convert matrix to vector and vice versa.
Download the test script test_hw5.m and the instructor's code yzSolve4X.p. You can run the test script without or with your code mySolve4X.m. First try on small n values (say, 100), then try on a large value (say, 1000). Submit your code, and output (on screen and graph, if any) from a small and a large n value. Your code should have a reasonable speed compared to the instructor's code. Write a paragraph to describe your observations.
(454, 30 pts) Problems NW (2nd ed.) 2.1, 2.2, 2.9.
(554, 30 pts) Problems NW (2nd ed.) 2.13 - 2.15.

HW 4: Due 2/7, Thur. (into instructor's mailbox, 1st floor mailroon, Dancun Hall)
(total points: 454/554 = 40/60) (Solution)

(10 pts) Program a Matlab function: function [x,iter] = myCGN(A,b,tol) that applies CG to the normal equations of Ax = b where A is nonsymmetric. Download the test script testCGN.m and the function mat4gmres.m. First run it for n = 1000 or less. Then run it for n = 5000 or larger. Submit your code and printout, and a short write-up (typed) explaining (not recounting) what you observe in this experiment.
(10 pts) Let A = aI + UV^T be non-singular and diagonalizable where U and V are m by r with r < m, and a is a scalar. For this kind of matrices A, prove that GMRES converges in at most r+1 iterations.
(10 pts) Let A = [e₂ e₃ ... e_m e₁] and b = e₁, where e_j is the j-th column of the identity (so A is a permutation of the identity). Prove that GMRES requires m iterations to solve Ax = b by looking at the Krylov subspaces generated by (A, b).
(554, 5 pts) Let A be defined as above. Derive the characteristic polynomial of A and then its eigenvalues.
(10 pts) Problems TB 38.6.
(554, 15pts) Exercises 5.3: problems 1-3 in the Krylov-Subspace Methods Notes.

HW 3: Due 1/31, Thur. (total points: 454/554 = 50/70) (Solution)

(20 pts) Implement the GMRES algorithm as is given in the instructor's "afternotes", with the Matlab interface function [x,iter] = myGMRES1(A,b,tol,maxit) where the output iter is the iteration number taken. Use the zero initial guess and the same stopping criterion as before. An adequate level of efficiency is always expected. It is a good practice to at first solve the inner least squares problem by a simple one-liner while make sure the Arnoldi part is done correctly. Then fully implement the whole algorithm.
Download the test script testGM1.m, the function fdm4pde1.m for matrix generation, and the instructor's GMRES code yzGMRES1.p. You can run the test script with or without your code. After you have written your code, start running the test script for a small m (say, m = 40) to debug and improve your code.
Submit your code and printout for runs with m = 80 or m = 100 (or even larger if you could handle it comfortably). Write a paragraph to summarize your observations.
(10 pts) Modify your (non-restarting) GMRES code into function [x,iter] = myGMRES1p(A,b,tol,maxit,M1,M2) so that it utilizes pre-conditioners M1 and M2 whenever they appear in the input list. The pre-conditioner is the product M = M1*M2 and the pre-conditioned system is inv(M)Ax = inv(M)b. See Matlab "help gmres" for more clues.
Download the test script testGM1p.m and the instructor's code yzGMRES1p.p. You can run the script without your code. When your code is ready, run the test script for n = 500 or larger (but use smaller values in the testing phase). Submit a copy of your program, the printput (graph not necessary), and a short write-up (typed) explaining what you observe in this experiment.
(10 pts) Modify your (non-restarting) GMRES code into a restarting GMRES code function [x,iter] = myGMRES2(A,b,p,tol,maxit) so that it restarts every p iterations. The input maxit is the maximum number of outer iterations allowed, and the output iter has 2 components, iter(1) is the number of outer iterations taken and iter(2) the number of inner ones taken in the last outer iteration. See Matlab "help gmres" for more clues.
Download the test scripts testGM2.m and testGM2a.m, and the instructor's code yzGMRES2.p. You can run these scripts without your code. When your code is ready, run testGM2 for m = 100 or larger (as large as you see fit), and also run testGM2a. Submit a copy of your program, the printputs (graph included) from your runs, and a short write-up (typed) explaining what you observe in the experiment (testGM2a).
(10 pts) Problems TB 35.5.
(554, 20 pts) (1) Problem TB 33.2. (2) Prove (10) in the "afternotes" on Krylov-subspace methods.

HW 2: Due 1/24, Thur. (total points: 454/554 = 30/50) (Solution)

(20 pts) Write a Matlab function: function [x,iter,info] = mySOR(A,b,omega,maxit,tol) for SOR method, where iter gives the iteration number taken, and omega is the parameter in SOR. Other input/output arguments have the same meanings as in HW 1. Please use the zero vector as your starting point and the same stopping criterion as in HW 1.
Download the test script test_SOR.m and the instructor's code zSOR.p. You can run test_SOR(1) or test_SOR(2), corresponding to 2 different matrices with default value at 1, with or without your code.
After you have debugged your code, run "test_SOR.m" with small N values and different ranges of omega to narrow down a good omega value (for example, you might try N = 50 and omega range = 0.8:0.1:1.9, and then 1.91:0.01:1.99 for matrix 1).
Run test_SOR(1) and test_SOR(2), using the best omega values that you have estimated, with much larger N values (as large as your code/computer could comfortably deal with).
Submit a copy of your program, a set of plots for initial runs, and the screen printout for your final runs. Write a paragraph to summarize your observations and thoughts.
(10 pts) Prove in detail that GS converges for any initial guess if the matrix A is symmetric positive definite.
(554, 20pts) Excercises 2 and 4 in the lecture notes on Stationary Iterative Methods.

HW 1: Due 1/17, Thur. (total points: 454/554 = 30/45) (Solution)

(20 pts) Write two Matlab functions:

function [x,info,res] = myJacobi(A,b,maxit,tol) function [x,info,res] = myGS(A,b,maxit,tol)
to implement the Jacobi and Gauss-Seidel methods, respectively, where maxit is the maximum number of iterations. Use the zero vector to start. The stopping criterion is that the residual 2-norm, norm(b-A*x), becomes less than tol*(1 + norm(b)). The output info is a character string indicating the status of termination: either "Converged at iteration ?" or "Maxit reached". And res is a vector containing residual norms, norm(b-A*x), for all iterations, counting iterations from 1.
Download the instructor's codes zJacobi.p, zGS.p and the test script test1.m. Read the test script.
Run test1.m for N = 50, 100, 200 (or larger if so desired), and with both method = 'Jacobi' and method = 'GS'. Submit your programs and output from these runs.
(10 pts) Prove in detail that the Jacobi method converges for any initial guess if the matrix A is column strictly diagonally dominant.
(554, 5 pts) Prove that the scheme x <-- x + α*(b-Ax) converges for any initial guess if the matrix A is symmetric positive definite and α>0 is less than 2 divided by the maximum eigenvalue of A.
(554, 10pts) Prove that in any matrix norm ||.||, ||A^k||^1/k ---> ρ(A) as k goes to infinity.