Numerical Computing Exercise 1: Inconsistent systems of equations [10 points] Solution

Consider the following inconsistent systems of equations:
(a) A1x = b1, with
1
(b) A2x = b2, with A2 = 0
1
1 1
1
2
0 0
1
1
1 2
b2 = 2 3 4
Find the least squares solution x∗ and compute the Euclidean norm of the residual, SE and RMSE.
Exercise 2: Polynomials models for least squares [20 points]
In this exercise, we consider two small datasets about the crude oil (crudeOil.txt) and kerosene (kerosene.txt) production by year in Europe in the period 1980-2012. By solving the following tasks, we will try to fit the data with different polynomial models and determine the best one.
(a) Write a Matlab function leastSquares.m which takes, as inputs, matrix A and vector b of a generic system Ax = b and returns, as outputs, the least squares solution x∗, the Euclidean norm of the residual, the SE and the RMSE. In order to test your function, write a script ex2a.m, in which you use leastSquares() to find the solutions of the two inconsistent systems of Exercise 1 and compare the results obtained with the ones you computed by hand.
(b) Consider the linear model yi = α1 + α2xi and apply it to the crude oil and kerosene production data in the period 1980-2011. Write a script linearModel.m in which you use leastSquares() to compute the least squares solution x∗ and the metrics of the residual. For each dataset, create a figure in which you plot the original data points and the linear model.
(c) Consider the quadratic model yi = α1 + α2xi + α3x2i and apply it to the crude oil and kerosene production data in the period 1980-2011. Write a script quadraticModel.m in which you use leastSquares() to compute the least squares solution x∗ and the metrics of the residual. For each dataset, create a figure in which you plot the original data points and the quadratic model.
(d) Consider the cubic model yi = α1 +α2xi +α3x2i +α4x3i and apply it to the crude oil and kerosene production data in the period 1980-2011. Write a script cubicModel.m in which you use leastSquares() to compute the least squares solution x∗ and the metrics of the residual. For each dataset, create a figure in which you plot the original data points and the cubic model.
1
(e) Compare the linear, quadratic and cubic models on the basis of the quality metrics computed above, by creating a table containing the results for the two models. Which one of the three models would you pick for the crude oil data? And for the kerosene? Provide an estimate of the crude oil and kerosene production in 2012 by using the three models and compare the values obtained with the real values reported in the data source. Comment on your results.
Exercise 3: Analysis of periodic data [20 points]
(b) Repeat the same analysis and plots of the previous point for both time series, by using the periodic model yi = α1 + α2 cos(2πxi) + α3 sin(2πxi) + α4 cos(4πxi) in the script periodicB.m.
(c) Compare the models of point (a) and (b). Was it beneficial to include more data? Which model would you prefer? Are you satisfied with the results obtained? If necessary, what would you suggest to improve your models? Motivate your answers.
Exercise 4: Data linearization and Levenberg-Marquardt method for the exponential model [20 points]
The file nuclear.txt contains the data on the nuclear electric power consumption by year in China in the period 1999-2006. We consider the power law model, expressed as:
.
(a) Find the least squares best fit by using data linearization and compute the RMSE both of the log-linearized model and of the original exponential model. Include in your report all the computations and the necessary steps, as explained in the slides of the tutorial.
(b) Write a function levenbergMarquardt() in which you implement the Levenberg-Marquardt algorithm for solving nonlinear least squares problems. Following again the slides of the tutorial, show how you can formulate the problem in order to solve it with Levenberg-Marquardt method and compute analytically all the necessary quantities. Finally, write a script ex4b.m in which you use the function levenbergMarquardt() to fit the data points and compute the RMSE.
(c) Compare the results obtained in points (a) and (b), by extending the script ex4b.m to produce a plot of the original data points together with the two models. Which model would you choose?
Exercise 5: Tikhonov regularization [15 points]
Let us consider an ill-posed problem Ax = b, with A ∈ Rn×m and b ∈ Rm. In order to solve it, we introduce a parameter α ∈R and we use linear least squares with Tikhonov regularization, given by:
. (1) x
(a) Derive analytically the optimal solution x∗ of Eq. (1). What is the purpose of the parameter α and how should we proceed to choose its value?
(b) We now consider the Hilbert matrix H ∈Rn×n, with entries defined as follows:
,
for every i,j = 1,...,n. Write a Matlab script illposedHilbert.m in which you generate H for n = 50,100,200,300,400,500,1000 and solve the problem Hx = b, for b = Hxexact and xexact generated through the function rand(n,1). To make your results reproducible, reinitialize the random number generator to its startup configuration by adding rng(’default’) at the beginning of your script. Produce also two figures in which you plot: (I) the condition number of H (use cond() in Matlab) against n; (II) the norm of the error kxexact −xk2 against n.
(c) We now focus our attention on the case n = 100. Write a Matlab script regularizedHilbert.m in which you estimate the regularized solution xreg according to Eq. (1) by using at least 10 different values of the parameter α (explain your choice of the values). To visualize the results, produce two figures in which you plot: (I) the norm of the error kxexact −xregk2 against the values of α; (II) kHx−bk2 against kxk2 for the different values of α. Comment your results.
Quality of the code and of the report [15 points]
The highest possible score of each project is 85 points and up to 15 additional points can be awarded based on the quality of your report and code (maximum possible grade: 100 points). Your report should be a coherent document, structured according to the template provided on iCorsi. If there are theoretical questions, provide a complete and detailed answer. All figures must have a caption and must be of sufficient quality (include them either as .eps or .pdf). If you made a particular choice in your implementation that might be out of the ordinary, clearly explain it in the report. The code you submit must be self-contained and executable, and must include the set-up for all the results that you obtained and listed in your report. It has to be readable and, if particularly complicated, well-commented.
Additional notes and submission details

Lecturer: Dr. E. Vecchi
In-class assistance
If you experience difficulties in solving the problems above, please join us in class either on Tuesdays 15:30-17:00 or on Wednesdays 13:30-15:00 in room C1.03.