Tutorial1D gaussian deconvolutionTo demonstrate how to use our package, we'll walk through the process of applying ADCG to 1D gaussian deconvolution. In this problem, the task is to decompose a 1d function into a sum of a few shifted gaussian functions. In other words, we get noisy observations of a function at a prespecified collection of points (often a 1d grid) and our task is to approximate as a sum of a few weighted gaussians: . In some applications, the function really is a sum of a few shifted gaussians, in which case we hope to recover their positions. Even this simple model has applications in multireturn lidar. In pictures, we want to take observations like this: Here the observed locations are a simple grid of points. We want to return estimates of and like this: Above we've plotted the estimated locations and weights as vertical lines on top of the denoised signal: . Because this is simulated data, we can plot the true locations to see how we did: 1d gaussian deconvolution as a sparse inverse problemRecall the definition of the sparse inverse problem: sparse inverse problem
We need to specify , , , , and . We'll denote the locations at which we observe as (in the plots above this would the 40 equally spaced points between 5 and 5) and the observed values as . Then we can take . The last object we have to specify is . . using SparseInverseProblems.jlAccording to the ADCG paper, we'll need to write functions to evaluate a few things: In terms of the julia package, this means we'll have to write an implementation of the ForwardOperator abstract type. That might be difficult. For problems with parameter spaces that are hypercubes in , that don't need to be extremely fast, we've provided a default subclass of ForwardOperator. If you have a complicated example, or you want a faster implementation, you might have to look at the ForwardOperator type. Our default abstract subtype is BoxConstrainedDifferentiableModel, for which we need to implement the following functions: from BoxConstrainedDifferentiableModel.jl
abstract BoxConstrainedDifferentiableModel <: ForwardModel #compute the measurement model psi(model :: BoxConstrainedDifferentiableModel, theta :: Vector{Float64}) #compute the jacobian of the forward model dpsi(model :: BoxConstrainedDifferentiableModel, theta :: Vector{Float64}) # Initial starting point for continuous optimization for the FW step. # Should return a good guess for $\arg\min_\theta \langle \psi(theta), v \rangle.$ # Often computed using a grid. getStartingPoint(model :: BoxConstrainedDifferentiableModel, v :: Vector{Float64}) = error("getStartingPoint not implemented for model $(typeof(model)).") # Box constraints on the parameters. # Returns a tuple of two vectors : lower bounds and upper bounds. parameterBounds(model :: BoxConstrainedDifferentiableModel) The only function here that isn't selfexplanatory is getStartingPoint. In order to solve the linear minimization problem: BoxConstrainedDifferentiableModel uses gradient descent — but as the function in this problem is nonconvex it needs a good initialization. This is what getStartingPoint provides. For simple examples with low dimensional parameter spaces, this is easy enough to do with gridding. implementationWe'll call our concrete subtype of BoxConstrainedDifferentiableModel SimpleExample for now. An instance of SimpleExample will need to know the points , so a first definition of SimpleExample might be: immutable SimpleExample <: BoxConstrainedDifferentiableModel evaluation_points :: Vector{Float64} end We'll first implement psi and dpsi: psi and dpsi
psf(theta, points) = exp((points . theta).^2/2.0) deriv_psf(theta, points) = exp((points . theta).^2/2.0).*(points . theta) psi(s :: SimpleExample, parameters :: Vector{Float64}) = psf(parameters,s.evaluation_points) dpsi(s :: SimpleExample, parameters :: Vector{Float64}) = reshape(deriv_psf(parameters,s.evaluation_points),length(s.evaluation_points),1) Not too hard. Next we'll provide an implementation of parameterBounds. Our single parameter is unconstrained: parameterBounds(model :: SimpleExample) = ([Inf], [Inf]) Finally we have to provide an initial starting point for the linear minimization step. The first step here is to modify the constructor for our type to take a grid of the parameter space and cache applied to those points: adding gridding to SimpleExample
immutable SimpleExample <: BoxConstrainedDifferentiableModel evaluation_points :: Vector{Float64} grid_points :: Vector{Float64} grid :: Matrix{Float64} SimpleExample(p,grid) = new(p,grid,psf(grid',p)) end Implementing getStartingPoint is then easy: we simply check all of our grid points and return the one that minimizes the objective: getStartingPoint(model :: SimpleExample, v :: Vector{Float64}) = [model.grid_points[indmin(model.grid'*v)]] That's everything. The final implementation is supplied with the package at ~.julia/v0.x/SparseInverseProblems/examples/simple_example , but we'll reproduce it here: simple_example.jl
using SparseInverseProblems import SparseInverseProblems: getStartingPoint, parameterBounds, psi, dpsi immutable SimpleExample <: BoxConstrainedDifferentiableModel evaluation_points :: Vector{Float64} grid_points :: Vector{Float64} grid :: Matrix{Float64} SimpleExample(p,grid) = new(p,grid,psf(grid',p)) end psf(theta, points) = exp((points . theta).^2/2.0) deriv_psf(theta, points) = exp((points . theta).^2/2.0).*(points . theta) psi(s :: SimpleExample, parameters :: Vector{Float64}) = psf(parameters,s.evaluation_points) dpsi(s :: SimpleExample, parameters :: Vector{Float64}) = reshape(deriv_psf(parameters,s.evaluation_points),length(s.evaluation_points),1) getStartingPoint(model :: SimpleExample, v :: Vector{Float64}) = [model.grid_points[indmin(model.grid'*v)]] parameterBounds(model :: SimpleExample) = ([Inf], [Inf]) a simple scriptThe final step is to simulate some data and finally call ADCG with our new ForwardModel. The following script generates some noisy data and applies ADCG to the output. It outputs three figures like those above. example script
include("simple_example.jl") using Gadfly, Compose evaluation_points = 5.0:0.25:5.0 model = SimpleExample(evaluation_points, 10.0:0.5:10.0) means = [2.0,1.0,3.5] k = length(means) weights = ones(k)+randn(k)*0.1 weights[1] = 0.2 target = max(phi(model,means',weights) + randn(length(evaluation_points))*0.1,0.0); function callback(old_thetas,thetas, weights,output,old_obj_val) #evalute current OV new_obj_val,t = loss(LSLoss(), output  target) println("gap = $(old_obj_val  new_obj_val)") return old_obj_val  new_obj_val < 1E1 end (means_est,weights_est) = ADCG(model, LSLoss(), target, Inf; callback=callback) #draw the observations alone draw(SVG("observations.svg", 5inch, 2inch), plot(x=evaluation_points,y=target)) #draw the observations with the locations and weights of the true blurs anno = Guide.annotation( compose(context(), circle([means,means], [zeros(k),weights], [2mm]), fill(nothing), stroke("blue"))) draw(SVG("truth.svg", 5inch, 2inch), plot(x=evaluation_points,y=target,anno)) #draw the estimation along with the predicted means/observations anno = Guide.annotation( compose(context(), circle([means_est',means_est'], [zeros(k),weights_est], [2mm]), fill(nothing), stroke("red"))) output_est = phi(model, means_est, weights_est) draw(SVG("est.svg", 5inch, 2inch), plot(x=evaluation_points, y=output_est,anno)) The only part of this script worth investigating in detail is the actual call to ADCG: (means_est,weights_est) = ADCG(model, LSLoss(), target, Inf; callback=callback) The returned values are pretty clear: a tuple of the estimated points along with their weights. The first three arugments tell ADCG which kind of problem we're solving, and what the loss function is. The last positional arugment gives the value of , which here we're setting to infinity. This tells ADCG that we want to approximately solve the nonnegative least squares problem. The mysterious part of this line is the optional argument callback. Recall that a sparse inverse problem requires a parameter that constrains the total mass of the returned solution, and in practice controls the number of sources we find. Often we don't know how to set this parameter, in which case we've found that an early stopping heuristic works very well: function callback(old_thetas,thetas, weights,output,old_obj_val) new_obj_val,t = loss(LSLoss(), output  target) println("gap = $(old_obj_val  new_obj_val)") return old_obj_val  new_obj_val < 1E1 end In each iteration of ADCG, we add a single source. The callback function simply tells ADCG to stop as soon as the drop in objective value from adding a new source is less than . To play around a bit more, try making the following modifications to the script and checking how they affect the output:
