Homework 4

Nonlinear least squares and source localization

\[ \def\argmin{\operatorname*{argmin}} \def\Ball{\mathbf{B}} \def\bmat#1{\begin{bmatrix}#1\end{bmatrix}} \def\Diag{\mathbf{Diag}} \def\half{\tfrac12} \def\ip#1{\langle #1 \rangle} \def\maxim{\mathop{\hbox{\rm maximize}}} \def\maximize#1{\displaystyle\maxim_{#1}} \def\minim{\mathop{\hbox{\rm minimize}}} \def\minimize#1{\displaystyle\minim_{#1}} \def\norm#1{\|#1\|} \def\Null{{\mathbf{null}}} \def\proj{\mathbf{proj}} \def\R{\mathbb R} \def\Rn{\R^n} \def\rank{\mathbf{rank}} \def\range{{\mathbf{range}}} \def\span{{\mathbf{span}}} \def\st{\hbox{\rm subject to}} \def\T{^\intercal} \def\textt#1{\quad\text{#1}\quad} \def\trace{\mathbf{trace}} \]

Homework submission instructions

• See the syllabus for the due date
• Follow the homework submission instructions

using LinearAlgebra, ForwardDiff, Plots, Random, Statistics
Random.seed!(1234); # seed for reproducibility

Consider the source localization problem

\[ \min_{x \in \Rn}\ \sum_{i=1}^m (\|x - c_i\|^2 - d_i^2 )^2, \tag{1}\]

where \(c_i\in\Rn\) represents the coordinates of the \(i\)th beacon and \(d_i\) is a positive number that represents the approximate distance from the unknown position \(x\in\R^2\) of the source to the \(i\)th sensor. Here we assume that \(n=2\), i.e., the source is localized in the plane, though this problem formulation extends to three or higher dimensions.

1 Gauss-Newton method

Implement the Gauss-Newton method for solving the smooth nonlinear least squares problem \[ \min_{x\in\Rn}\ \half\|r(x)\|^2, \quad r:\Rn\to\R^m. \tag{2}\] The resulting function should have the following signature:

"""
Gauss-Newton method for nonlinear least-squares using a
backtracking linesearch (Armijo rule).

Inputs:
  `r`: residual function (function from `x` to `R^m`)
  `J`: Jacobian of `r` (function from `x` to `R^(m,n)`
  `x0`: initial guess
  `maxits`: maximum number of iterations
  `ϵ`: stopping tolerance
  `μ`: backtracking parameter (default 1e-4)`

Outputs:
  `x`: approximate solution
  `k`: number of iterations
  `err`: vector of gradient norms at each iteration

Algorithm stops when the 2-norm of the gradient of the
nonlinear least-squares objective 0.5∥r(x)∥² is less
than `ϵ` or the number of iterations exceeds `maxits`.
"""
function gauss_newton(r, J, x0; maxits=100, ϵ=1e-6, μ=1e-4)
    # ...
    return x, k, err
end

2 Solve the source localization problem

This function generates an instance of the source localization problem, where m beacons are placed at random, and the source is placed at a random location. The function also generates distances—corrupted by noise—from each beacon to the unknown source location.

Data generator for source localization problem

"""
Structure for the source localization problem with `m` beacons on the plane.
"""
struct SourceLocalization
    c::Matrix{Float64} # 2x`m` matrix with the beacon positions
    d::Vector{Float64} # vector of distances (noisy)
    x::Vector{Float64} # true position (you're not meant to know this!)
end;

"""
Generate an instance of `SourceLocalization` that contains data
for `m` beacon locations, target location, and distance data
between the beceans and target.

Inputs:
  `m`: number of beacons
  `η`: noise level (default 0.1)

Output:
    `SourceLocalization` object
"""
function generate_data(m; η=0.1)
    c = randn(2, m)  # beacon positions
    x = randn(2)     # true (unkown) position
    d = [norm(x - ci) + η*randn() for ci in eachcol(c)]
    return SourceLocalization(c, d, x)
end;

1. Write a function localize_source that solves a Gauss-Newton method to solve the source localization problem described by Equation 1 and use it to solve the following tasks:

"""
Localize source using Gauss-Newton method
Input:
  `sl`: SourceLocalization object

Outputs:
    `x`: estimated position
    `itn_error`: vector of gradient norms at each iteration
"""
function localize_source(sl::SourceLocalization)
    # code here ...
    return x, itn_error
end;

2. Solve an instance of the source localization problem for m = 10 beacons with a noise level of \(\eta=0.1\). That is, execute these commands:

sl = generate_data(10)
x, _ = localize_source(sl);

3. On a single scatter plot, visualize the positions of the beacons, your predicted guess deduced from the beacon range data, and the position of the true (unknown) position. (In practice, of course, we couldn’t visualize the true position, but it’s useful here for our experiments.)

You’re welcome to use the code below to plot the positions of the beacon, the true (unknown) position, and a predicted guess deduced from the beacon data.

Position plotter

function plotmap(sl::SourceLocalization, xguess)
  scatter(leg=:outertopright,frame=:box, aspect_ratio=:equal)
  scatter!(sl.c[1,:], sl.c[2,:], label="beacons", shape=:square, ms=7)
  scatter!(sl.x[1,:], sl.x[2,:], label="true position", shape=:xcross, ms=7)
  scatter!(xguess[1,:], xguess[2,:], label="guess", c="orange", ms=7)
end;

3 Solution confidence

Because we’re placing the beacons randomly, we cannot precisely predict the quality of the source position that we compute from a particular configuration of beacon positions. We might intuit, however, that the accuracy of our estimated location improves with additional beacons. That raises the question: how does the error in the estimated position depend on the number of beacons? In partial answer to this, write code that carries out the following experiment:

• Localize the source for 200 trials for a fixed value of \(m\).
• Determine the mean and standard deviation for those 200 random trials.
• Repeat the trials for \(m=10, 20, 30, 40, 50, 60\) beacons.
• Plot the mean and standard deviation of these trials as a function of \(m\).

4 Gauss-Newton vs gradient descent

For \(m=10\) beacons, plot the per-iteration progress in the norm of the gradient of the nonlinear least-squares objective Equation 2 for the Gauss-Newton method and for the gradient descent method with Armijo backtracking linesearch.

5 Beck, Problem 5.5(i)

(Beck, second edition)