COMPENG 3SK3
Project 2
1 Background
Newton’s Method in Optimization
Course ID: Instructor: Due date:
COMPENG 3SK3 February 26, 2023
To meet safety requirements transportation authorities need to monitor the conditions of highway and railway tunnels and detect any structural changes. One can use LiDAR (light detection and ranging) to measure and model the tunnel interior surface. For this purpose, laser reflection markers are painted on the tunnel surface so that LiDAR can measure their 3D positions. The challenge is that the LiDAR readings are noisy, which affects the precision of the 3D measurements. One way to suppress the sensor noises is to scan the markers with LiDAR for K > 1 times, and combine the data in an optimal fashion to improve the measurement precision.
In this project, you are to apply Newton’s method to perform the above optimization task. A detailed derivation of the algorithm is presented in the next section.
2 Method
2.1 Notations
pi = (xi , yi , zi )T is the ground-truth position of the i-th marker on the tunnel surface.
pei is the estimate of pi made by the proposed algorithm.
qk = (ak , bk , ck )T is the position of the LiDAR detector in the k-th observation.
T is the noisy measurement made in the k-th observation. ik ikikik
dik = ||p?ik ?qk||2 = p(p?ik ? qk)T (p?ik ? qk) is the measured distance between the i-th marker and the LiDAR position in the k-th observation.
N is the number of markers.
K is the number of observations.
2.2 Illustation
Figure 1 illustrates how our LiDAR-based tunnel monitoring system works. There are N markers {pi}Ni=1 on the
tunnel surface. The LiDAR observation positions {qk}Kk=1 are evenly spaced on a straight line. These positional data
are stored in a K × 3 matrix (given to you in ’observation R5 L40 N100 K21.mat’). In the k-th observation, the
LiDAR makes noisy measurements of the N markers, denoted by {p?ik}Ni=1. The assemble of these K measurements
are recorded in a K × N × 3 data array (given to you in ’pts R5 L40 N100 K21.mat’). Besides, the distances
{dik}N,K between the i-th marker and the k-th LiDAR observation point are stored in an N × K matrix (given i =1,k =1
to you in ’dist R5 L40 N100 K21.mat’).
For each marker, K noisy coordinates and K distances, one for each of the K observations, are available to
estimate the position of this marker. Simply averaging K noisy coordinates can help reducing the noise to some extent, but this is not optimal. In this project, Newton’s method will be applied to obtain an optimal estimate of marker. In next subsection, how to solve this problem with Newton’s method is discussed in detail.
2.3 Formulation
For simplicity, we drop the index of marker i from the subscript of pi, pei, p?ik and dik in the following discussions.
2023 Winter Term
Figure 1: A demonstration for the system prototype.
Given the measured distances {dk}Kk=1 between {qk}Kk=1 and the marker p and the noisy coordinates {p?k}Kk=1, the error can be represented as
where λ is a constant to control the importance of the two error terms.
The optimal coordinate of the marker is 1KK
p2
Let r ∈ RK in Eq.(3) denote the vector of residuals.(1) can be rewritten into
The Jacobian matrix of r with respect to p can be represented as
where the k-th row of J is, its gradient with respect to p can be obtained by the rule [1]
The first item in Eq.(7) is equivalent to
Then, Eq.(7) can be rewritten into a more concise form
Next, according to the chain rule in calculus, the gradient of the loss function Eq.(4) with respect to p is
g = JT r + λ X 2(p ? p?k).
k=1
Regarding the Hessian matrix of loss function Eq.(4), JT J can approximate the Hessian matrix of 1 PK
2 k=1k and the Hessian of ||p k||2 is 3 × 3 identity matrix I3×3. Then, the Hessian of loss function Eq.(4) is
H ≈ JT J + 2λKI3×3 (10) With g and H, we can obtain the descent direction for Newton’s method. (11) In each iteration, the coordinate will be updated according to
p(j+1) = p(j) ? αH?1g (12)
where α is used to control the step size along the descent direction. α can be either calculated by line search algorithm [3] or set to a constant. Line search is beyond the scope of this course, so simply set α to a constant (e.g., α = 0.1).
3 Work to do
Your task is to implement Newton’s method by writing your own program according to the above formulation. This Newton’s method based denoising algorithm is point-wise, i.e., it takes the corresponding measurements of one marker as input and predicts the optimal coordinate of this marker. So you should process each marker one at a time using the same algorithm1, and put these predictions together into an N × 3 matrix.
Requirements:
Implement the above-mentioned Newton’s method based denoising algorithm with Matlab. Don’t use any built-in functions such as fminunc and cvx. Please use only transpose, norm, addition, subtraction, multiplication, and division (including inv). Python is also acceptable. If you choose Python, please use basic operators of numpy to write your own programs and do not use any advanced library functions such as scipy.optimize and pytorch.
Submit your code with detailed documentation.
1Use the same λ, initialization scheme p(0) and termination condition to process each marker.
Use pseudo code to describe all necessary details of your implementation.
Evaluate the quality of estimated coordinates of markers by computing the root mean square error (RMSE) [4]
between your result {pei}Ni=1 and the ground-truth markers {pi}Ni=1
vu u 1 XN
100
Fine-tune the initialization p(0)
– Let λ = λop and keep λ fixed, sketch three curves of RMSE against iterations given the following
initialization schemes:
average the K measured coordinates;
choose one of the measured coordinates;
set it as a random vector.
– Comment on the needed number of iterations to converge for the above three cases.
Individual work. Comments:
The ground-truth coordinates of markers are organized as a N × 3 matrix in ’gt R5 L40 N100 K21.mat’. The given ground-truth markers can only be used to compute the RMSE of your result {pei}Ni=1 for the purpose of self-evaluation.
Your code will be verified on another dataset which is unknown to students.
A template for your implementation is shown as follows4
RMSE = tN
– Find λop that can make the converged RMSE as small as possible2;
Fine-tune λ
– Plot the curve of converged RMSE against λ. The range of λ is [ 1 λop,100λop]3.
Algorithm 1: Newton’s method based denoising algorithm input : {p?ik }N,K , {di k }N,K , {qk }K
Extract measured coordinates and distances corresponding to the i-th marker; Initialization;
while not converged do
update the coordinate of the i-th marker according to the Newton’s method; Record pei;
References
[1] Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. Technical University of Denmark, 7(15):510, 2008.
[2] Wikipedia contributors. Gauss–newton algorithm — Wikipedia, the free encyclopedia, 2023. [Online; accessed 2-February-2023].
2Use the average of K measured coordinates (i.e., K1 PKk=1 ?pk) as the initialization p(0). Make sure that the converged RMSE is
smaller than the RMSE obtained by simply average K1 PKk=1 ?pk, and the reference for RMSE reduction is 0.0068.
3λ?=0and0<λop <1.
4The while loop can be replaced by a for loop (for it=1:max iterations) where max iterations is the maximum number of iterations
and defined by you. If the termination condition for Newton’s method is satisfied, break this for loop.