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. • pˆ = (xˆ , yˆ , zˆ )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 1KK X(||p−qk||2 −dk)2 +λX||p−pˆk||2 (1) k=1 k=1 E(p)= 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. r1 ||p−q1||2 −d1  r2 ||p−q2||2 −d2  2 pe=argmin X(||p−qk||2 −dk)2 +λX||p−pˆk||2. (2) k=1 k=1 p(p−q1)T(p−q1)−d1  p(p−q2)T(p−q2)−d2  r= . = . = . . (3) .. .  rK ||p−qK||2 −dK p(p−qK)T(p−qK)−dK Then, Eq.(1) can be rewritten into 1KK Xrk2 +λX||p−pˆk||2. (4) k=1 k=1 E(p)= The Jacobian matrix of r with respect to p can be represented as 2 J=∂x ∂y ∂z (5)  . . .   . . .  ∂rK ∂rK ∂rK ∂x ∂y ∂z ∂r1 ∂r1 ∂r1  ∂x ∂y ∂z ∂r2 ∂r2 ∂r2  2 where the k-th row of J is h∂rk ∂rk ∂rk i (p−qk)T hx−ak y−bk z−ck i ∂x ∂y ∂z = p(p−qk)T(p−qk) = rk+dk rk+dk rK+dk As for ||p − pˆk||2, its gradient with respect to p can be obtained by the rule [1] . (6) (7) (8) (9) r 2 [2] g = The first item in Eq.(7) is equivalent to r ∂rk + λ k ∂y  ∂rk k=1 ∂z 2(p − pˆ ). k Then, Eq.(7) can be rewritten into a more concise form k=1 ∂r1 ∂r2 ∂x ∂x ∂r1 ∂r2 ∂y ∂y ∂r1 ∂r2 ∂z ∂z r  ··· ∂rK 1 ∂||x||2 = 2x. ∂x Next, according to the chain rule in calculus, the gradient of the loss function Eq.(4) with respect to p is K ∂rk  K X∂x X ∂x r2 ··· ∂rK .=JTr ∂y . ··· ∂rK . ∂z rK K 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 − 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 ∆p = −H−1g. (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. 3 • 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 i=1 ||pei − pi||2. (13) 1 2 3 4 5 6 i=1,k=1 i=1,k=1 k=1 output: {pei}Ni=1 fori←1toNdo 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. 4 [3] Wikipedia contributors. Line search — Wikipedia, the free encyclopedia, 2022. [Online; accessed 2-February- 2023]. [4] Wikipedia contributors. Root-mean-square deviation — Wikipedia, the free encyclopedia, 2023. [Online; accessed 2-February-2023].