Sparse regularization has proven to be an indispensable tool in many areas, including inverse problems  and compressive sensing [2, 3]. If a signal is known to have a sparse or compressible (quickly decaying) representation , this information can be used to formulate optimization problems of the form
where is a measurement matrix used to measure the true signal ,
is a vector of data, andis a threshold that depends on the characteristics of measurement error. In compressive sensing, it is possible to obtain recovery guarantees given properties of the true signal and . In more general inverse problems, these guarantees have not been found; it is therefore appropriate to consider (1) as a regularization approach to the least squares problem. For example, in the seismic setting, where (1) has been particularly useful , is a linearized Born-scattering operator, is the curvelet transform, and is seismic data. While there are several popular algorithms that solve (1), e.g. SPArsa , the SPG  algorithm has been particularly useful for seismic imaging [7, 8, 9].
Many inverse problems contain unknown nuisance parameters that must be estimated in order to recover the solution. In seismic imaging, the source wavelet is typically unknown. The main contribution of this paper is to extend the approach of  to the sparse inversion context, and derive simple modifications of standard sparse solvers to incorporate solutions of unknown nuisance parameters on the fly.
The paper proceeds as follows. In section 2, we introduce the seismic imaging problem with unknown wavelet, and formulate it as an extended sparsity promoting optimization problem (3). In section 3, we review the ideas recently proposed in  that allow nuisance parameters to be estimated on the fly, and show how to incorporate these ideas into existing sparsity promoting formulations. In section 4, we develop an extended SPG algorithm to solve (3), and we present numerical results in section 5.
2 Imaging with Unknown Wavelet
Seismic imaging is an approach to obtain a gridded subsurface velocity perturbation from seismic data, given a smooth starting model. Experiments are conducted by placing explosive sources on the surface and recording the reflected waves with an array of receivers on the surface. The data,
, in this case represents the Fourier transform of the recorded time series for frequency. The corresponding modeling operator, , defines a linear relation between the recorded data for the frequency and the velocity perturbation. The statistical model for data given is
where is a statistical model for the measurement error, which is typically modeled as Gaussian, and are unknown complex source wavelet coefficients. Note that the model (2) is no longer linear in the decision variables —it is bilinear. Since the perturbation is known to be sparse in the Curvelet frame , formulation (1) has been successfully used to recover  when the source wavelet is known. In full generality, the joint inverse problem for the perturbation and wavelet is given by
3 Variable Projection
We begin by considering the problem
If we define , problem (4) is of the form
where for any given , one can easily find
using the convenient formula
This is basically a generalization of the variable projection algorithm .
4 Projected Regularized Inversion
In the previous section, we showed how to solve the extended problem (4). However, the formulation (3) is more important to us from the modeling perspective, since it is always easier to provide a noise threshold than to figure out the ‘right’ sparsity level . In fact, the SPG algorithm solves formulation (1) by solving a series of subproblems that find the sparsity level automatically given an input threshold .
Suppose we find such that . Can we expect that the corresponding minimizers of (4) coincide with the minimizers of (3)? This question is answered in surprising generality by [12, Theorem 2.1]: as long as any minimizer of (3) satisfies , then the set of minimizers of and match, and where .
This general result points to using the following strategy: solve by Newton’s method
Evaluating is straightforward: we simply use the projected gradient method detailed in section 3. However, is more difficult, since the most general variational results for value functions  require linearity of the forward model, which is violated by (2). Nonetheless, if we treat as fixed, then by [12, Theorem 6.2], we get
where solves (4) for .
We note that the expression above is an approximation to the derivative, and the quality of the approximation remains to be determined. If the source weight can be estimated fairly quickly (so that it is not changing significantly between iterations), the approximation above becomes exact. For the experiments in the next section, we found that the proposed Newton iteration gives nearly the same result as the one with a fixed, ‘true’ source-weight. We also verified that when we pick that is reachable within our computational budget of 150 iterations, the algorithm correctly finds the root ; see figure 4 (b).
5 Numerical Results
For the experiments we use a Matlab framework for seismic imaging and modelling , and the CurveLab toolbox . Both of these are freely available for non-commercial purposes. The algorithm to solve 3 is based on the SPG code , which is also available for download.
We generate data for the velocity perturbation defined on a 201 x 301 grid with 10 m spacing depicted in figure 1 for 6 frequencies between 5 and 25 Hz, 301 equispaced receivers and 15 composite sources, all located at the top of the model. We note that his leads to a underdetermined problem with 27090 equations and 60501 unknowns. We use SPG to solve (3) either with fixed or with estimated using the procedure outlined above. Since there is no noise in this example we use and run the algorithm for a fixed number of iterations (100 in this case).
Note that there is a fundamental non-uniqueness in the problem; if we multiply the source-weights with a constant factor, we can compensate for this by dividing the reconstructed model by the same factor. Therefore, we normalize the results such that the source-weights for each reconstruction have the same norm (i.e., is the same for all reconstructions).
The reference result using the true source signature is shown in figure 2 (a). If we do not estimate the source signature and use , we do not get a good reconstruction, as is shown in figure 2 (b). Finally, if we estimate the source signature according to the strategy outlined in this paper, we obtain the result depicted in figure 2 (c).
6 Discussion and Conclusions
We have proposed a novel method for estimating nuisance parameters in the context of sparsity regularized inverse problems, and in particular we have focused on source wavelet estimation in seismic imaging. The method draws on the idea of variable projection in order to estimate nuisance parameters on the fly, and can be implemented via a straightforward modification to existing sparse solvers.
Numerical experiments demonstrate that the source wavelet can be recovered successfully in this manner (figure 3), and that the recovery of primary parameters (specifically of the image) is improved when the wavelet is estimated (figure 2).
Note that after section 3, we can already solve (4), but nonetheless lot of effort is devoted in section 4 to develop a method for solving (3). The main point here is that while it is difficult to come up with a reasonable value for in (4), it is straightforward to come up with a good value for in (3). In fact, given a finite computational budget, and no estimate for , one can always pick and perform a fixed number of iterations. In this mode, the algorithm in section 4 solves several (4) problems inexactly, picking the corresponding sequence of values according to iteration (11). This is exactly what was done to obtain the numerical examples in section 5.
-  Felix J. Herrmann, Peyman P. Moghaddam, and Chris C. Stolk, “Sparsity- and continuity-promoting seismic imaging with curvelet frames,” Journal of Applied and Computational Harmonic Analysis, vol. 24, no. 2, pp. 150–173, 2008, doi:10.1016/j.acha.2007.06.007.
-  D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
-  E. J. Candès, L. Demanet, D. L. Donoho, and L. Ying, “Fast discrete curvelet transforms,” Multiscale Modeling and Simulation, vol. 5, pp. 861–899, 2006.
-  Aleksandr Y. Aravkin, Xiang Li, and Felix J. Herrmann, “Fast seismic imaging for marine data,” in ICASSP. ICASSP, 2012, ICASSP.
-  S.J. Wright, R.D. Nowak, and M.A.T. Figueiredo, “Sparse reconstruction by separable approximation,” Signal Processing, IEEE Transactions on, vol. 57, no. 7, pp. 2479 –2493, july 2009.
-  E. van den Berg and M. P. Friedlander, “Probing the pareto frontier for basis pursuit solutions,” SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 890–912, 2008.
Felix Herrmann, Michael Friedlander, and Ozgur Yilmaz,
“Fighting the Curse of Dimensionality: Compressive Sensing in Exploration Seismology,”IEEE Signal Processing Magazine, vol. 29, no. 3, pp. 88–100, May 2012.
-  Felix J. Herrmann and Xiang Li, “Efficient least-squares imaging with sparsity promotion and compressive sensing,” Geophysical Prospecting, vol. 60, no. 4, pp. 696–712, 2012.
-  Xiang Li, Aleksandr Y. Aravkin, Tristan van Leeuwen, and Felix J. Herrmann, “Fast randomized full-waveform inversion with compressive sensing,” Geophysics, vol. 77, no. 3, pp. A13–A17, 2012.
-  Aleksandr Y Aravkin and Tristan van Leeuwen, “Estimating nuisance parameters in inverse problems,” Inverse Problems, vol. 28, no. 11, pp. 115016, 2012.
-  G Golub and V Pereyra, “Separable nonlinear least squares: the variable projection method and its applications,” Inverse Problems, vol. 19, no. 2, pp. R1, 2003.
-  Aleksandr Y. Aravkin, James V. Burke, and Michael P. Friedlander, “Variational properties of value functions,” submitted to SIAM J. Optimization, arXiv:1211.3724, 2012.
-  T. van Leeuwen, “A parallel matrix-free framework for frequency-domain seismic modelling, imaging and inversion in matlab,” 2012, https://www.slim.eos.ubc.ca/releases.
-  David Donoho Lexing Ying Emmanuel Candes, Laurent Demanet, “CurveLab,” 2006, http://www.curvelet.org/.
-  E. van den Berg and M. P. Friedlander, “SPGL1: A solver for large-scale sparse reconstruction,” June 2007, http://www.cs.ubc.ca/labs/scl/spgl1.