Development of a spectral source inverse model by using
generalized polynomial chaos
Kyongmin Yeoa,∗, Youngdeok Hwangb , Xiao Liuc , Jayant Kalagnanama
a

arXiv:1801.03009v1 [physics.ao-ph] 9 Jan 2018

IBM T.J. Watson Research Center, Yorktown Heights, NY 10598, USA
b
Department of Statistics, Sungkyunkwan University, Seoul, Korea
c
Department of Industrial Engineering, University of Arkansas, Fayetteville, AR, USA

Abstract
We present a spectral inverse model to estimate a smooth source function
from a limited number of observations for an advection-diffusion problem.
A standard least-square inverse model is formulated by using a set of Gaussian radial basis functions (GRBF) on a rectangular mesh system. Here, the
choice of the collocation points is modeled as a random variable and the generalized polynomial chaos (gPC) expansion is used to represent the random
mesh system. It is shown that the convolution of gPC and GRBF provides
hierarchical basis functions for a spectral source estimation. We propose a
mixed l1 and l2 regularization to exploit the hierarchical nature of the basis polynomials to find a sparse solution. The spectral inverse model has
an advantage over the standard least-square inverse model when the number of data is limited. It is shown that the spectral inverse model provides
a good estimate of the source function even when the number of unknown
parameters (m) is much larger the number of data (n), e.g., m/n > 50.
Keywords: Inverse model, Advection-diffusion equation, Source estimation,
Stochastic programming, Generalized polynomial chaos, Uncertainty
quantification
1. Introduction
Air pollution, generated by either anthropogenic or natural causes, poses
a major public health threat. Not only long-term (Hoek et al., 2013), but also
acute exposure (Phalen and Phalen, 2011) over a certain threshold can cause
health problems. Due to its immense importance, there have been substantial development in the computational modeling of the transport of air-borne
∗

kyongmin@gmail.com

Preprint submitted to Elsevier

January 10, 2018

pollution over the past decade (Byun and Schere, 2006; El-Harbawl, 2013;
Fast et al., 2006). However, prediction of air pollution by using these computational models requires extensive prior information on the distribution and
magnitudes of pollution emission sources, which in most cases is incomplete
or has high uncertainty (Thunis et al., 2016). Moreover, in many cases, it is
of greater interest to identify the source of pollution when abnormally high
pollution is observed in the air quality monitoring network to mitigate a
possible public health hazard. This atmospheric inverse problem to find the
pollution emission source using a set of measurements from a sensor network,
has attracted significant attention in the atmospheric science community.
One of the fundamental building blocks of the inverse model is the atmospheric dispersion process, which is modeled by an advection-diffusion
equation (Stockie, 2011). Deterministic approaches, in which a partialdifferential-equation constrained optimization problem is solved to minimize
a convex loss function, typically a square difference between the observations
and the computational model predictions, have been adopted from the field
of atmospheric data assimilation and used widely (Eckhardt et al., 2008; Issartel et al., 2007; Martinez-Camara et al., 2014; Pudykiewicz, 1998). Since
the optimization formulation for an inverse model usually leads to underdetermined or ill-conditioned system, much research effort is focused on regularizing the solution. Recently, inverse models exploiting Bayesian inference
have become popular (Chow et al., 2008; Keats et al., 2007; Rajaona et al.,
2015), due to the strength of the Bayesian methods in dealing with noisy
and incomplete data. However, Bayesian inference usually requires extensive number of simulations to compute the posterior distribution, such as
the Markov-chain Monte Carlo method. In the pollution source inversion
problem, by exploiting the linearity of the advection-diffusion equation, the
adjoint advection-diffusion operator can be used to reduce the computational
cost(Keats et al., 2007). More general approaches to mitigate the high computational cost have been proposed by either accelerating the convergence
of a Monte Carlo simulation (Marzouk et al., 2007), constructing surrogate
models (Li and Marzouk, 2014), or developing a low-dimensional representation (Lieberman et al., 2010; Roosta-Khorasani et al., 2014).
Most of the previous inverse models consider either estimating the magnitudes of the source at each computational grid points by combining a large
volume of heterogenous data (de Foy et al., 2015; Hwang et al., 2018; Issartel
et al., 2007), or finding the locations and magnitudes of one or a few point
sources from a limited number of data (Keats et al., 2007; Marzouk et al.,
2007). In this paper, we propose an inverse model based on a regularized
optimization formulation to estimate a smooth source function from a small
number of observations. First, we follow the conventional approach of ap2

proximating a smooth function by a set of Gaussian radial basis functions
centered at a set of collocation points of a rectangular mesh system. Obviously, the solution of the inverse model is strongly dependent on the choice
of the collocation points. To relax such dependency on the mesh system, we
propose modeling the choice of the collocation set as a random variable. A
stochastic inverse model is formulated on this random mesh system and the
generalized Polynomial Chaos expansion (gPC) (Xiu, 2007) is employed to
tackle the stochastic inverse problem. It is shown that the stochastic formulation results in a spectral inverse model, in which the unknown smooth
function is approximated by hierarchical polynomials. The spectral inverse
model has an advantage over the conventional least-square inverse model,
particularly when the number of data is limited, due to its capability of
hp-refinement (Karniadakis and Sherwin, 2005).
This paper is organized as follows. Section 2 describes a least-square formulation of the advection-diffusion problem by using an adjoint operator. In
section 3, we reformulate the deterministic least-square problem as a stochastic problem by using gPC. In section 4, a mixed l1 - and l2 -regularization is
introduced to exploit the hierarchical nature of the basis polynomials and
an algorithm to solve the optimization problem is presented. The proposed
inverse model is tested in section 5. Finally, the concluding remarks are given
in section 6.
2. Least-square inverse model
2.1. Forward model
We consider the following advection-diffusion problem,


Aφ(x, t) = Q(x) x ∈ D
n · ∇φ(x, t) = 0 x ∈ ∂Dout


φ(x, t) = 0
x ∈ ∂Din

(1)

in which D is a rectangular domain D = (0, L1 ) × (0, L2 ) with the boundary
∂D and n is an outward normal vector on ∂D = ∂Dout ∪ ∂Din . The outflow
and inflow boundaries are ∂Dout = {x : x ∈ ∂D, n · u(x, t) ≥ 0} and
∂Din = {x : x ∈ ∂D, n · u(x, t) < 0}. The advection-diffusion operator is
defined as


∂
∂
∂
∂
(uj (x, t)φ(x, t)) −
Kij (x, t)
φ(x, t) ,
Aφ(x, t) = φ(x, t) +
∂t
∂xj
∂xi
∂xj
(2)
in which u(x, t) is the fluid velocity, and K(x, t) is a symmetric secondorder tensor for the diffusivity. We assume that u(x, t) and K(x, t) are
3

smooth functions with uniformly bounded derivatives of all orders. The
source strength Q(x) is the unknown function but assumed to be smooth.
Furthermore, we consider a non-negative source, i.e., Q(x) ≥ 0 for every
x ∈ D.
Contrary to the usual computational prediction problem, where equation
(1) is solved for known Q(x) to estimate φ(x, t) at the sensor locations, the
inverse model aims to estimate Q(x) from the given observations Φ(to ) at
the observation time to . Here, instead of estimating Q(x) directly by solving
an infinite dimensional optimization problem, Q(x) is approximated by the
sum of a set of basis functions to reduce the dimensionality of the problem.
Let Q∗ (x) be the finite-dimensional approximation,
∗

Q (x) =

Nk
X

βi Pi (x).

(3)

i=1

Here, P is a set of basis functions, Nk is the total number of the basis
functions and β = (β1 , . . . , βNk ) denotes the coefficients. The problem of
estimating a continuous surface is reduced to a problem of finding Nk coefficients, β. There are many possible choices forP
the basis functions as long as
Q∗ (x) satisfies the non-negativity condition:
Pi (x)βi ≥ 0. Here, a set of
Gaussian radial basis functions (GRBF) located at the collocation points of
a rectangular mesh is used as the basis functions;


1 |x − y i |2
1
exp −
,
(4)
Pi (x) =
2π(c∆)2
2 (c∆)2
in which ∆ is the distance between the neighboring collocation points, y i is
the location of i-th collocation point and c is an O(0.1) parameter. Obviously,
in the limit of Nk → ∞, P converges to the Dirac delta function Pi (x) →
δ(x − y i ) and, hence, Q∗ (x) converges uniformly to Q(x). Using GRBF, the
non-negativity condition becomes, simply, βi ≥ 0 for i = 1, · · · , Nk .
Now, the advection-diffusion equation can be written as
Aφ(x, t) =

Nk
X

βi Pi (x).

(5)

i=1

Since A is a linear operator, we can exploit the superposition of the solutions.
Let φi (x, t; βi ) be the solution for the i-th GRBF;
Aφi (x, t; βi ) = βi Pi (x).
Then, clearly,
φ(x, t) =

Nk
X
i=1

4

φi (x, t; βi ).

Moreover, from the linearity of A,
φi (x, t; βi ) = φei (x, t)βi .

(6)

Here, φei is the solution for Pi with a unit strength, i.e., βi = 1.
The computational model output is related to the observation by an inner
product with respect to a sensor function (χ) as
Φi (to ) = hφ(x, t), χ(x, t; xoi , to )i + i for i = 1, · · · , No ,

(7)

in which Φi (to ) is the observation at the i-th sensor, xoi is the location of
the i-th sensor, to is the time of the measurement, No is the total number
of the sensors, and  is a Gaussian white noise representing the errors in the
measurement as well as the computational model. The angle bracket denotes
an inner product
Z
Z 0
ha(x, t), b(x, t)i =
a(x, t + τ )b(x, t + τ )dτ dx.
x∈D

−∞

The sensor function depends on the types of the sensor or the data used in
the analysis. In this study, the sensor function is defined as
χ(x, t; xoi , to ) =

1
δ(x − xoi ){H(t − to + Tχ ) − H(t − to )}.
2Tχ

(8)

Here H(t) is a Heaviside function, which is zero for t < 0 and one otherwise,
and Tχ is an time-average window of the sensor.
By comparing (6) and (7), β can be related to the measurement by
f + .
Φ = Xβ

(9)

Here, ΦT = (Φ1 (to ), · · · , ΦNo (to ))T and
eij = h χ(x, t; xo , to ), φej (x, t) i.
X
i

(10)

f ∈ RNo ×Nk relates β to the observation Φ. Since
The dispersion matrix X
f can be computed by solving exactly the same partial differential equation
X
(2) for each GRBF (Pi ), the same numerical solver can be recycled.
The source strength Q∗ (x) can be obtained by finding β from the following least-square minimization problem
1
f 2 + R(β),
arg min kΦ − Xβk
2
2
N
β∈R k ,β≥0

(11)

f is
in which R(β) is a regularization. Since we consider the case Nk  No , X
a rank-deficient matrix and a regularization is required to enforce uniqueness
of the solution.
5

2.2. Adjoint model
f requires to solve the advectionIt is important to note that computing X
diffusion equation for Nk times, which makes it computationally impractical
as Nk becomes large. Moreover, as will be discussed in section 3, when
the model uncertainty is considered, the total number of computation easily
blows up to O(103 ∼ 104 ). To circumvent these difficulties, an adjoint model
is employed in this study. Reducing the number of repetitive computations
from the number of GRBF, Nk , to the number of observations, No , an adjoint
model is computationally more tractable (Keats et al., 2007).
Here, the adjoint model is briefly described. Define a conjugate field (φ∗ )
as
(
hφ(x, t), χ(x, t; xoi , to )i = hφ∗i (x, t; to ), Q(x)i for i = 1, · · · , No
(12)
for t ≥ to
φ∗i (x, t) = 0
Then, the adjoint operator is obtained from the Lagrangian duality relation;
hAφ, φ∗i i = hφ, A∗ φ∗i i = hφ(x, t), χ(x, t; xoi , to )i,

(13)

which gives
A∗ φ∗i (x, t; to ) = −

∂ ∗
φ −u·∇φ∗i −∇·(K ·∇φ∗i ) = χi , for t ∈ (−∞, to ), (14)
∂t i

for χi = χ(x, t; xoi , to ). The adjoint model (14) is solved backward in time
from to . For more details, see Pudykiewicz (1998).
Once the i-th conjugate field φ∗i is computed by solving (14) with appropriate boundary conditions (Hourdin et al., 2006), φ at the i-th sensor is
computed as
hφ(x, t), χ(x, t; xoi , to )i = hφ∗i (x, t; to ), Q∗ (x)i =

Nk
X

hφ∗i (x, t; to ), Pj (x)iβj .

j=1

(15)
Repeating the process for No conjugate fields, the observation vector is
Φ = Xβ + ,

(16)

in which
Xij = hφ∗i (x, t; to ), Pj (x)i.
f = X. The coefficients β are computed by solving
It is trivial to show that X
the same least-square minimization problem (11).
6

Figure 1: Example of the collocation points (black dot) and a basis function Pi (x).

3. Generalized polynomial chaos for model uncertainty
The Gaussian radial basis function, P, distributes the source strength
β computed from (11) in the space centered on the collocation points of
a rectangular mesh system. Define a rectangular mesh system by a tensor
product of two one-dimensional collocation sets,
W = s1 ⊗ s2 .

(17)

Here, s1 and s2 denote the sets of collocation points in the x1 - and x2 directions, respectively;
s1 = {xi1 ; xi1 = x01 + i × ∆, i = 1, · · · , Nk1 },
s2 = {xi2 ; xi2 = x02 + i × ∆, i = 1, · · · , Nk2 },
in which x01 and x02 are, respectively, the left and bottom ends of the mesh
system, and Nkj is the number of the collocation points in the j-direction,
Nk = Nk1 × Nk2 . An example of W and GRBF is shown in figure 1.
The solution Q∗ (x) of the standard least-square (LS) inverse model depends on the choice of W. For example, if a local peak of Q(x) does not
coincide with one of the collocation points, the LS inverse model will result in a poor accuracy. In general, there is no standard rule of choosing
W. In this study, we propose to represent the uncertainty in the choice of
W as a random variable. Let W ∗ (ω) be a random variable with a uniform
distribution;
W ∗ (ω) = (s1 + ξ1 (ω)∆) ⊗ (s2 + ξ2 (ω)∆).
(18)
Here, ξ1 (ω) and ξ2 (ω) are real random variables defined over a probability
space (Ω, S, P), in which Ω is the sample space, S is the σ-algebra, P is the
7

probability measure, and ω is an element of Ω. Both ξ1 (ω) and ξ2 (ω) are
defined over Γ1 = Γ2 = (−0.5, 0.5) with the probability density functions,
ρi (ξ) : Γi → R+ , ρ1 (ξ1 ) = ρ2 (ξ2 ) = 1, i.e. ξi ∼ U(−0.5, 0.5) for i = 1, 2. Note
that W ∗ (ω) corresponds to a random translation of W, which uniformly
covers the entire computational domain.
On the random collocation system, W ∗ (ω), equation (16) becomes
Φ = X(ω)β(ω) + ,

(19)

in which
Xij (ω) = hφ∗i (x, t; to ), Pj (x; ω)i,
and



1 |x − (y i + ∆ξ(ω))|2
1
.
exp −
Pi (x; ω) =
2π(c∆)2
2
(c∆)2
Taking an expectation over ξ, (19) is
Φ = Eξ [X(ω)β(ω)] + ,
and a least-square inverse model can be formulated as
1
arg min kΦ − Eξ [X(ω)β(ω)]k22 + R(β(ω)), s.t. Eξ [Q∗ (x; ω)] ≥ 0 ∀x ∈ D,
2
β(ω)
(20)
in which
∗

Q (x; ω) =

Nk
X

βi (ω)Pi (x; ω).

i=1

Hereafter, the obvious dependence on ξ is omitted in the expectation, i.e.,
E[f ] = Eξ [f ].
Following Xiu (2007); Xiu and Karniadakis (2002), the gPC expansion is
employed to approximate the stochastic functions;
β(ξ) =

X(ξ) =

M
X
i=0
M
X

βbi Ψi (ξ),

(21)

ci Ψi (ξ).
X

(22)

i=0
i

Here, Ψ (ξ) is an orthonormal polynomial basis in a bivariate polynomial
space (Z2P ) constructed by a tensor product of one-dimensional polynomial
spaces
2
O
P
Z2 ≡
W i,P ,
i=1

8

bi (x1 ), for c = 0.1.
Figure 2: First four modes of the one-dimensional gPC coefficients, P

in which

W i,P ≡ v : Γi → R : v ∈ span{ψm (ξi )}Pm=0 .
As ξ is a uniform random variable, the Legendre polynomial is chosen as the
basis polynomial ψm (ξ) (Xiu and Karniadakis, 2002).
From the gPC approximation,
Z
M X
M
M
X
X
i
j
c βb
ci βbi , (23)
E[X(ξ)β(ξ)] =
X
Ψi (ξ)Ψj (ξ)ρ(ξ)dξ =
X
i=0 j=0

ξ∈Γ

i=0

in which Γ = Γ1 × Γ2 , ρ(ξ) = ρ1 (ξ1 ) = ρ2 (ξ2 ) = 1, and M denotes the total
number of the basis functions excluding the mean (zero-th order) component,
M = (P + 1)2 − 1. Similarly,
!
Nk
M
X
X
bj (x) ,
(24)
E[Q∗ (x; ξ)] =
βbj P
i

i=1

i

j=0

and
i
i
bm
blm
X
= hφ∗l (x, t; to ), P
(x)i.

bi (x), is
Here, the gPC coefficient for GRBF, P


Z
1 |x − (y j + ∆ξ)|2
1
i
b
Pj (x) =
exp −
Ψi (ξ)dξ.
2π(c∆)2
2
(c∆)2
9

(25)

(26)

b can be easily computed. For example, the first mode of
The coefficients P
the expansion is


Z
j
2
1
|x
−
(y
+
∆ξ)|
1
0
bj (x) =
Ψ0 (ξ)dξ
(27)
P
exp −
2
2
2π(c∆)
2
(c∆)
"
!
!#
2
yij − xi + 0.5∆
1 Y
yij − xi − 0.5∆
√
√
erf
=
− erf
.
4∆2 i=1
2c
2c
b
bi (x). It is shown that P(x)
Figure 2 shows the first four modes of P
constitutes spatial hierarchical basis functions to approximate Q(x). The advantage of using the hierarchical basis functions over increasing the number of
GRBF (Nk ) is discussed in section 4.
4. Regularized optimization formulation
Using the gPC expansion, the minimization problem (20) becomes,
Nk X
M
M
X
X
1
i bi 2
b
c
bj (x) ≥ 0 ∀x ∈ D.
X β k2 + R(β), s.t.
βbij P
arg min kΦ −
i
2
N
(M
+1)
b
k
β∈R
i=0
i=1 j=0

(28)
b =
In the minimization problem, the number of parameters to estimate is dim(β)
Nk × (M + 1). Note that, if we use the forward simulation approach shown in
section 2.1, even for Nk = O(100) and P = 10, the total number of numerc becomes O(104 ). On the other hand, using
ical simulations to compute X
the adjoint model, the total number of the numerical simulations remains as
c can be computed efficiently by evaluating the
No , which is ∼ O(10) and X
b
inner product hφ∗ , P(x)i
with a numerical integration, such as a Gaussian
quadrature.
In equation (28), βb should satisfy the non-negativity constraint, E[Q∗ (x)] ≥
0 for every x ∈ D. Because it is difficult to directly impose the non-negativity
condition for every x, we propose an indirect constraint based on the stochastic collocation approximation (Xiu and Hesthaven, 2005). In the stochastic
collocation method, the stochastic functions are approximated by the Lagrange polynomials as
β(ξ) =
P(x; ξ) =

N
X
i=1
N
X
i=1

10

βei Li (ξ),

(29)

e i (x)Li (ξ).
P

(30)

Here, Li (ξ) is the i-th Lagrange polynomial, which is Li (ξj ) = δij for the
j−th collocation points in Γ, ξj . Again, L(ξ) is constructed by the tensor
product of two one-dimensional Lagrange polynomials in Γ1 and Γ2 . The
stochastic collocation method corresponds to a deterministic sampling and
the coefficients are easily computed by the function evaluations at each cole i (x) = P(x; ξi ). Since GRBF is
location points in Γ, i.e., βei = β(ξi ) and P
positivie, P(x; ξ) ≥ 0 for all (x, ξ) ∈ D × Γ, the non-negativity condition
implies β(ξi ) ≥ 0 for i = 1, · · · , N . In other words, in the stochastic collocation approach, it is sufficient to impose the non-negativity constraint only
on the parameters, not on the field. The stochastic collocation coefficients,
b as
β(ξi ), are related to the modal gPC coefficients, β,
β(ξi ) '

M
X

βbj Ψj (ξi ).

j=0

Then, the constraint on the smooth function surface E[Q∗ (x)] ≥ 0 for every
x ∈ D can be approximated by the following linear constraint
βi (ξj ) =

M
X

βbik Ψk (ξj ) ≥ 0 for i = 1, · · · , Nk & j = 1, · · · , N.

(31)

k=0

Or,
Lβb ≥ 0,

(32)

in which L is a block diagonal matrix for Ψi (ξj ).
As an analogy to the finite element analysis (Karniadakis and Sherwin,
2005), either h- or p-type refinement can be used to increase the resolution of
the approximation, Q∗ (x). Let Nk be the number of GRBFs for a reference
case. In the h-type refinement, the grid space is decreased as ∆q = ∆/q for
q ∈ N+ , which makes the total number of unknown parameters m = Nk × q 2 .
In the p-type refinement, Nk is fixed and the maximum order of Ψ(ξ) is
increased, which results in m = Nk × p2 for p = P + 1. Because of this
quadratic dependence, in both h- and p-type refinements, the number of unknown parameters can easily overwhelm the number of observations upon
a refinement. For example, in the numerical experiments in section 5, the
b is ∼ O(103 − 104 ), while the number of
number of unknown parameters, β,
data ∼ O(10). As a result, the optimization formulation results in a highly
ill-posed system, of which solution heavily relies on the choice of the regularization. To alleviate the difficulty, we develop a regularization strategy,
b
which exploits the hierarchical nature of the GRBF-gPC coefficients, P(x),

11

b 0 (x), repreAs shown in equation (27), the zero-th gPC mode of GRBF, P
j
sents the average source strength in (−0.5∆, 0.5∆) × (−0.5∆, 0.5∆) centered
b 0 (x) ≥ 0, the non-negativity constraint implies that
at y i . Because P
j
βbi0 ≥ 0 for i = 1, · · · , Nk .
At the same time, when the mean source strength, βbi0 , is zero, the fluctuations
around the mean, represented by the higher-order gPC modes, should also
be zero, i.e.,
βbik = 0 for every k ≥ 1, if βbi0 = 0.
Therefore, in the p-type refinement, the number of unknown parameters can
be effectively reduced by identifying non-zero elements in βb0 .
From these observations, we propose the following mixed l1 - and l2 regularizations,
b = λ1 kβb0 k1 + λ2 kβb0 k2 ,
R(β)
(33)
2
in which βb0 = ({βb1 }T , · · · , {βbM }T )T . The Least Absolute Shrinkage and
Selection Operator (LASSO), or l1 -regularization, is one of the most widely
used regularization methods to find such a “sparse” solution for the so-called
“large m, small n” problem (large number of parameters and small number
of data (Tibshirani, 1996). As discussed above, for the p-type refinement,
applying LASSO only for βb0 is enough to guarantee a sparse solution. Hence,
LASSO is applied only to the zeroth mode, βb0 , and the higher-order terms
are regularized by the standard Tikhonov regularization.
k
From the non-negativity constraint, βb0 ∈ RN
≥0 . Then, the regularized
optimization problem can be written as
1
arg min kΦ − E[X(ξ)β(ξ)]k22 + λ1 E[β(ξ)] + λ2 tr (Cov [β(ξ)]) ,
2
β(ξ)

(34)

s.t. E[Q∗ (x)] ≥ 0, ∀x ∈ D.
Here, it is clear that the first tuning parameter λ1 controls the sparsity in the
solution, β, and the second one λ2 prevents the overfitting by regularizing
the total variation around the mean.
Furthermore, to consider the spatial smoothness of Q(x), the fused LASSO
is used. In the fused LASSO, l1 -norm is imposed on the difference between a
directly connected parameters (Tibshirani et al., 2005). For example, in the
x1 -direction, the fused LASSO regularization is
X
kG1 βb0 k1 =
|βbi0 − βbj0 |,
(35)
(i,j)∈N1

12

in which N1 is an index set for βb0 directly connected in the x1 -direction, i.e.,
N1 = {(i, j) : (y i − y j ) · e1 = ∆ for i = 1, · · · , Nk − 1 and i < j ≤ Nk }.
In other words, the fused LASSO is equivalent to imposing l1 -penalty in the
gradient of Q∗ (x). We can define the difference matrix in the x2 -direction,
G2 , similar to G1 . Then, the regularization can be written as a mixed generalized LASSO (Tibshirani. and Taylor, 2011) and Tikhonov regularization,
b = λ1 kS βb0 k1 + λ2 kβb0 k2 ,
R(β)
2
in which

(36)




γI
S = G1  .
G2

The coefficient γ determines the relative weight of the standard LASSO to
the fused LASSO.
Finally, the optimization problem for the spectral source reconstruction
is
M

X
1
ci βbi k2 + λ1 kS βb0 k1 + λ2 kβb0 k2 , s.t. Lβb ≥ 0.
X
arg min kΦ −
2
2
2
N
(M
+1)
b
k
β∈R
i=0

(37)

This optimization problem is solved by using the Alternating Direction Method
of Multipliers (ADMM). ADMM is one of the most widely used method for
a large-scale optimization (Boyd et al., 2010).
Rewriting (37) as a constraint optimization problem,
1
b 2 + λ1 kαk1 + λ2 kβb0 k2 , s.t. F βb = ζ,
kΦ − X βk
2
2
2
N
(M
+1)
b
k
β∈R
arg min

(38)

c0 , · · · , X
cM ], F = [{S, 0}> , L> ]> , ζ = (α> , θ > )> , and θ =
where X = [X
b + . Note that F is padded with a null matrix {0}, because the general(Lβ)
ized LASSO, S, is applied only to βb0 . The subscript + indicates a projection
onto a positive set R≥0 , i.e., (b)+ = max(b, 0). The augmented Lagrangian
form is
b 2 + λ1 kαk1 + λ2 kβb0 k2
b α, v) = 1 kΦ − X βk
L(β,
2
2
2
ω
+ ωv > (F βb − ζ) + kF βb − ζk22 ,
2

(39)

in which ω is a constant. Then, the solution procedure for the minimization
problem is as follows;
13

1. Set initial conditions for βb(0) and other variables.
e −1 (X > Φ)}+ .
βb(0) = {(X > X + ωF > F + λ2 I)
b whose first Nk
Here, Ie is an identity matrix of the dimension of β,
elements are zero, which imposes the Tikhonov regularization on βb0 .
Then,
α(0) = S βb0, (0) ,
θ (0) = (Lβb(0) )+ ,
v (0) = 0.
2. For t > 0, update βb(t+1) by
e −1 (X > Φ + ωF > (ζ (t) − v (t) )).
βb(t+1) = (X > X + ωF > F + λ2 I)
3. Next, update ζ, i.e., α and θ by
α(t+1) = sign([F βb(t+1) + v (t) ]U ) × (1 − (λ1 /ω)/k[F βb(t+1) + v (t) ]U k2 )+ ,
θ (t+1) = ([F βb(t+1) + v (t) ]L )+ ,
 (t+1) 
α
(t+1)
ζ
=
.
θ (t+1)
Here, v is a vector of dimension Nk + N , i.e., the number of GRBFs
and the number of the quadrature points for the evaluation of the
non-negativity condition. The operator [·]U takes only the upper Nk
elements of the vector for the generalized LASSO and [·]L subsets the
lower N elements of the vector for the non-negativity constraint.
4. Finally, update v,
v (t+1) = v (t) + (F βb(t+1) − ζ (t+1) ).
5. Repeat steps 2 ∼ 4, until the improvement
kβb(t+1) − βb(t) k2 < 0
for a pre-specified tolerance level 0 .
The Tikhonov regularization is imposed in Step 2. In Step 3, the fused
LASSO is imposed on α and the non-negativity constraint is imposed by a
projection on the non-negative set in θ.

14

Figure 3: (a) The computational domain
and the location of the sensors. (b) Examples
R
of a few time-integrated adjoint field, φ∗i (x, t; t0 )dt.

5. Numerical experiments
We use numerical simulations to study the behavior of the spectral inverse model and the results are compared with the standard least-square
inverse models with LASSO and fused-LASSO regularizations (Tibshirani
et al., 2005; Tibshirani, 1996). The computational domain for this case study
is D = (−10, 35) × (−10, 35). The adjoint equation (14) is numerically integrated for T = (−5, 0) by using a third-order low-storage Runge-Kutta
method and an upwind finite volume method. The computational grid and
time step sizes are chosen as δx1 = δx2 = 0.2 and δt = 1/600.
Figure 3 (a) shows the computational domain D as well as the locations
of the sensors. There are total 36 sensors in D. The sensors are located at
a rectangular mesh of which grid size is ∆χ = 5. The time averaging window of the sensor in equation (8) is Tχ = 1 and the time of observation is
t0 = 0. The computational parameters are chosen similar to real operational
air pollution measurement conditions. In an operational air pollution measurement, usually a high frequency sensor measurement is averaged over one
to three hours to remove measurement noise.
A few adjoint functions are shown in figure 3 (b). To mimic the atmospheric dispersion process, the wind field is generated by a Fourier transform
of a set of Ornstein-Uhlenbeck processes. The velocity in the x1 -direction is
!
h X
h
X
u1 (x, t) = U10 (t) + Real
(b
ur (l, m, t) + ib
ui (l, m, t))eik(l,m)·x , (40)
l=1 m=1

in which h is the maximum number of the Fourier modes and k(l, m) =
( L2π1 l, L2π2 m) is the wavenumber. The Fourier coefficients are computed by
15

Figure 4: (a) True source surface together with the sensor numbers. (b) The observation
(Φ) at each sensors.

solving the following Langevin equation;
u
br,i (l, m, t)
δt + S(l, m)δW,
(41)
TL
in which TL is a relaxation timescale, S is a scale parameter, and δW
q denotes
the Wiener process ∼ N (0, δt). In this study, TL = 2 and S = 2 TL (l22+m2 )
δb
ur,i (l, m, t) = −

are used (Pope, 2000). The velocity in the x2 -direction is computed from the
divergence-free condition;
∂u1 (x, t) ∂u2 (x, t)
+
= 0.
∂x1
∂x2

(42)

In other words,
!
p
p
X
X
L2 l
u2 (x, t) = U20 (t) − Real
(b
ur (l, m, t) + ib
ui (l, m, t))eik(l,m)·x .
L1 m
l=1 m=1
(43)
0
(t) are also obtained by solving the same Langevin
The mean components U1,2
q
equation, but with S = 5 T2L . The diffusivity is computed by an isotropic
Smagorinsky model, which is typically used in the atmospheric dynamics
models (Byun and Schere, 2006):
Kij (x, t) = Kh (x, t)δij ,
s
 

∂u1 (x, t) ∂u2 (x, t)
∂u1 (x, t) ∂u2 (x, t)
2
Kh (x, t) = (Cs ∆s )
−
+
+
.
∂x1
∂x2
∂x2
∂x1
Here, δij is the Kronecker delta, Cs (= 0.1) is the Smagorinsky coefficient,
and the length scale ∆s = max(L1 , L2 )/2πh.
16

Figure 5: Estimated source surface by (a) gPC-LASSO, (b) fused LASSO, and (c) LASSO.
The black dots indicate the centers of GRBF and the white crosses in (a) are the locations
of the sensors.

5.1. Case study 1
Figure 4 (a) shows the true source surface, Q(x). The data is a snapshot
of the concentration observation from the sensor network, shown in figure
4 (b). The goal is to estimate the source surface Q(x) from the snapshot
observation.
Figure 5 shows the estimated source surface, Q∗ (x), by the spectral inverse model (gPC-LASSO) and the LS inverse model with fused LASSO, and
LASSO regularizations. Hereafter, we use fused LASSO and LASSO to refer
the LS inverse model with fused LASSO and LASSO regularizations, respectively. In gPC-LASSO, the grid space of the collocation set W is ∆ = 5 and
the scale parameter of GRBF, P(x), is set to c = 0.25. The maximum order
of the Legendre polynomial is P = 5 in each direction. For fused LASSO
and LASSO, ∆ = 2 and c = 0.5 are used. The penalty parameter for l1
regularization is set to λ1 = 10−2 for the all three models. For gPC-LASSO,
the penalty parameter for l2 regularization is λ2 = 10−6 . In gPC-LASSO and
fused LASSO, γ = 0.5 is used.The ratio of the number of unknown parameters to the data is NkgP C × p2 /36 ' 52 for gPC-LASSO and Nk /36 ' 13 for
the LS inverse model. So, gPC-LASSO has many more unknown parameters
than the standard LS inverse model. We aimed to choose the computational
parameters of the gPC-LASSO and the LS inverse models, such as ∆ and P ,
to make the ratio of the unknown parameters to the data size roughly equal
for a fair comparison. However, it is found that the standard LS inverse
model showed a very poor performance when the number of parameters is
increased by the h-type refinement. So, we do not show the LS inverse model
results for ∆ = 1, where Nk /36 ' 52.
In figure 5 (a), it is shown that gPC-LASSO provides a very good approximation to the true source surface. For a quantitative comparison, we
17

b
Figure 6: The modal coefficients β.

define a normalized l2 error,
R
eQ =

{Q∗ (x) − Q(x)}2 dx
R
.
Q2 (x)dx

(44)

The normalized l2 errors are 0.03, 0.17, and 0.47 for gPC-LASSO, fused
LASSO, and LASSO, respectively.
In gPC-LASSO, although the grid space of W is bigger than the size of
the true emission source, the estimated source surface is very close to the true
surface. Due to the hierarchical nature of the basis functions, gPC-LASSO
can provide a good estimate of the function surface even when the local
maximum of Q(x) is located in the middle of the collocation points. While
the fused LASSO correctly identifies the spatial pattern and the magnitude of
the source (5 b), Q∗ (x) depends strongly on the choice of the collocation set
W. LASSO identifies the locations of large Q(x), but fails to provide a good
approximation of the function surface. The ratio of maximum value of the
source surface, {max Q∗ (x)}/{max Q(x)}, is 0.99 and 0.93 for gPC-LASSO
and fused LASSO, respectively, while that of LASSO is 1.46, as LASSO tries
to fit the data with a fewer number of stronger sources.
b Among the total
Figure 6 shows the values of the modal coefficient β.
1,764 parameters, there are only about 280 non-zero coefficients. It is shown
that the magnitude of βb decreases as the order of polynomial increases, i.e.
going right on the horizontal axis. The magnitudes of βb of the highest order
polynomials are very small, indicating that the maximum order of Legendre
polynomial used in this case study is enough to resolve Q(x).
In figure 7, gPC-LASSO is tested for two different collocation sets. From
W used in figure 5 (a), in figure 7 (a), W is shifted in the x1 -direction by 0.5∆
18

Figure 7: The source surfaces from gPC-LASSO for two different collocation sets.

Figure 8: The source surfaces from gPC-LASSO with the grid space ∆ = 8 and the
maximum order of Legendre polynomial of (a) P = 4 and (b) P = 8.

and, in figure 7 (b), by 0.5∆ in the x2 -direction. As expected, it is shown
that gPC-LASSO is not sensitive to the choice of the basis collocation set.
The errors for the two new collocation sets are roughly the same, eQ ' 0.04,
which is very similar to the error of the reference case (figure 5 a), eQ ' 0.03.
Figure 8 shows Q∗ (x) from gPC-LASSO at a coarser grid resolution,
∆ = 8. The same penalty parameters are used λ1 = 10−2 and λ2 = 10−6 ,
while the scale parameter is changed, c = 0.15. Figure 8 (a, b) shows the
estimated source surfaces for two different maximum order of gPC expansion;
P = 4 and 8. For P = 4, gPC-LASSO under-resolves the source surface.
When P is increased to 8 (fig. 8 b), even with the large grid space, gPCLASSO is able to approximate Q(x) fairly well. The l2 error changes from
eQ = 0.26 for P = 4 to eQ = 0.09 for P = 8. It is also observed that
increasing P further does not increase the accuracy of the source estimation.
19

Figure 9: Changes of eQ with respect to the penalty parameters; (a) λ1 and (b) λ2 . In
(a), the solid symbol (•) denotes ∆ = 5 and the hollow symbol (◦) is for ∆ = 8.

In gPC-LASSO, there are two penalty parameters; λ1 and λ2 . As discussed in section 4, the l1 penalty parameter λ1 controls the sparsity of
the mean E[β], while λ2 is related with the spectral decay of the variance,
tr(Cov(β)). To show the effects of these penalty parameters on the solution,
eQ is computed for a wide range of λ1 ad λ2 . In figure 9 (a), eQ is shown as
a function of λ1 for ∆ = 5 and 8 for a fixed λ2 = 10−6 . In general, eQ is not
very sensitive to λ1 as long as λ1 is sufficiently small < 10−2 . For ∆ = 5,
eQ seems to have a local minimum around λ1 = 10−2 , which then increases
rapidly for a larger λ1 . It is worthwhile to note that decreasing λ1 does not
have a significant effect on the solution, implying that the LASSO regularization does not play an important role in imposing sparsity in the solution.
It is well known that, in a least-square regression problem, non-negativity
constraint alone, without LASSO, is enough to recover sparsity (Buckstein
et al., 2008; Slawski and Hein, 2011; Wang et al., 2011). In gPC-LASSO,
although non-negativity constraint is not directly imposed on βb0 , the linear
constraint (Lβb ≥ 0) effectively imposes the non-negativity constraint on the
mean components, βb0 , which explains why the generalized LASSO regularization for the mean component, kS βb0 k1 , does not have a significant effect
on the solution.
The effect of λ2 is shown in figure 9 (b). Again, eQ is not sensitive to λ2 .
However, when λ2 becomes larger than 10−3 , the higher-order modes of βb
are significantly suppressed, which results in the rapid increase of eQ .
Another important parameter of gPC-LASSO is the lengthscale of the
basis Guassian envelope, c∆. The scale parameter, c, essentially controls the
smoothness of the solution. Figure 10 shows the effects of c on the square

20

Figure 10: Effects of the scale parameter, c, on Q∗ (x). The hollow symbols (◦) denote the
square loss function (left axis) and the solid symbols (•) are for eQ (right axis).

loss function,
V (c) = kΦ −

M
X

ci (c)βbi k2 ,
X
2

i=0

together with the changes in eQ . It is shown that eQ has a local minimum
around c = 0.25 , while V (c) is almost constant up to c = 0.25, which starts
to increase rapidly afterward. As c is changed from 0.25 to 0.28, V (c) is
increased by more than two orders of magnitude. For small c, gPC-LASSO
tries to approximate Q(x) with sharp, peaked polynomials, resulting in a
highly oscillatory surface. Hence, although gPC-LASSO is able to find a
solution to faithfully fit Φ, eQ becomes large due to the high oscillation. On
the other hand, for larger c, the basis polynomials become too smooth to
approximate Q(x), which makes both V (c) and eQ grow. Figure 10 provides
a hint to select the scale parameter without a priori knowledge by finding
the maximum c which satisfies V (c) < δc for a threshold level δc .
In figure 11, gPC-LASSO is tested against noisy observation data. The
36 observations shown in figure 4 (b) are perturbed by a multiplicative noise;
Φ∗ = Φ × (1 + 0.05%),

(45)

in which % is a uniform random variable in (−1, 1). Due to the small number
of data (36 observations), an inverse model can be very sensitive to the noise.
However, it is observed that, even with the noise in the data, gPC-LASSO
can satisfactorily identify the spatial pattern of the true source surface Q(x).
Unlike the noise-free case, Q∗ (x)s from two different grid space (∆ = 5 and 8)
have very similar accuracy; eQ = 0.24 and 0.25 for ∆ = 5 and 8, respectively.
We compute the ensemble error
Z
E[eQ ] = eQ (%)ρ(%)d%.
21

Figure 11: The estimated source surfaces from noisy data; (a) ∆ = 5 and (b) ∆ = 8.

The integral is computed by a Latin hypercube sampling (McKay et al., 1979)
with 200 samples. This sampling method is known to achieve better accuracy
in numerical integration. The ensemble error for ∆ = 8 is E[eQ ] = 0.27. We
performed the same test with LASSO, which is found to be very sensitive to
the noise due to the large parameter-to-data ratio.
5.2. Case study 2
In the second case study, the emission surface, Q(x), has a larger lengthscale than the grid space ∆ (figure 12 a), and the emission surface extends
to the area not covered by the sensors. Figure 12 (b) shows Q∗ (x) from
gPC-LASSO. The parameters are chosen the same with the first case study,
∆ = 5, P = 5, λ1 = 10−2 , and λ2 = 10−6 , except for the scaling parameter
c = 0.5. It is shown that gPC-LASSO again provides a good approximation
of Q(x). The normalized l2 -error is eQ = 0.04 for gPC-LASSO. In figure 12
(c) and (d), Q∗ (x) from LASSO and fused LASSO are shown, respectively.
The model parameters are ∆ = 5, c = 0.5, λ1 = 10−2 , and λ2 = 10−6 . Note
that the model parameters for fused LASSO and LASSO are changed from
the first case study, because using a finer mesh, ∆ = 3 as in the first case
study, results in a large error. The normalized l2 errors are 0.13 and 0.14,
respectively, for LASSO and fused LASSO.
Figure 13 shows the behavior of eQ with respect to the penalty parameters, λ1 and λ2 . Similar to the previous results (figure 9), it is shown that eQ
is not sensitive to λ1 and λ2 . Although the normalized l2 error, eQ , shows a
local minimum around λ2 = 10−1 , the difference between the local minimum
and eQ at smaller λ2 is only 0.01.
Figure 14 shows E[eQ ] for noisy data. Again, the 36 observations, Φ,
are perturbed by the same multiplicative noise as in equation (45) and the
22

Figure 12: (a) True source surface and the estimated source surfaces by (b) gPC-LASSO,
(c) LASSO, and (d) fused LASSO. The white crosses in (a) indicate the locations of the
sensors.

Figure 13: Dependence of eQ on the penalty parameters λ1 (•) and λ2 (◦).

23

Figure 14: The expectation of eQ as a function of λ1 for noisy data. The error bar denotes
one standard deviation.

expectation is computed by the Latin hypercube sampling with 200 samples.
Contrary to the noise-free case, the l1 regularization has a large impact on the
accuracy of Q∗ (x). When λ1 is very small, gPC-LASSO finds a solution which
overfits the noisy data. In other words, gPC-LASSO results in a solution with
a highly oscillatory surface to minimize the square loss function, kΦ − Xβk22 ,
for the noisy Φ. For a larger λ1 , the smoothness of the solution is recovered by
the generalized LASSO regularization, which makes eQ smaller. It is shown
that, as λ1 increases, the standard deviation of eQ is decreased, implying that
the generalized LASSO regularization is robust to the noise.
6. Summary
In this paper, we present a spectral inverse model to estimate a source
function from a limited number of data for an advection-diffusion problem.
One of the conventional methods of approximating a smooth source function,
Q(x), is to discretize the computational domain by a set of collocation points
(W) and compute the coefficients, β(W), of a basis function, such as GRBF.
However, in such a collocation approach, the estimated function surface,
Q∗ (x), strongly depends on the choice of W. To remove the dependence on
the fixed collocation set, we formulate a stochastic least-square inverse model
on a random mesh system, W ∗ (ω). The generalized polynomial chaos expansion (gPC) is employed to approximate the resulting stochastic functions;
the source-receptor relation X(ω) and the source strength β(ω). By using
gPC, a spectral inverse model is formulated, where Q(x) is approximated
by hierarchical polynomials. The spectral approach has advantages over the
conventional collocation-based method in that Q∗ (x) is not as strongly dependent on W, and the spatial sparsity in Q(x) is more effectively recovered.
The non-negativity constraint of Q(x) ≥ 0, ∀x ∈ D, is replaced by a linear
constraint, Lβb ≥ 0, by comparing gPC with a stochastic collocation method.
24

Finally, a mixed l1 and l2 regularization is proposed based on the hierarchical
nature of the basis polynomials. An ADMM algorithm is presented to solve
the regularized optimization problem.
The solution behavior of the proposed method (gPC-LASSO) is investigated for two case studies and the model error is compared with the collocation methods with two most widely used regularization methods, LASSO
and fused LASSO. It is confirmed that gPC-LASSO is not sensitive to the
choice of W and provides a very good approximation to the source surface
even when the number of unknown parameters is more than 50 times larger
than the number of data. It is shown that gPC-LASSO outperforms both
LASSO and fused LASSO. For the noise-free data, the regularization does
not play an important role for gPC-LASSO because the non-negativity constraint alone is enough to explore the sparsity in the solution. However,
when noisy is added to the observations, the regularization provides a more
accurate approximation by enforcing the smoothness in the solution.
In this paper, we consider a two-dimensional steady source. However,
the same framework can be easily extended to a more general situations of
three-dimensional and / or unsteady source. It is worthwhile to note that,
although our method is presented in the context of air pollution problem,
the framework can be applied to a general function estimation problem for a
smooth, non-negative function.
References
Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimization
and statistical learning via the alternating direction method of multipliers.
Foundations and Trends in Machine Learning 2010;3(1):1–122.
Buckstein A, Elad M, Zibulevsky M. On the uniqueness of nonnegative
sparse solutions to underdetermined systems of equations. IEEE Trans Inf
Theory 2008;54:4813–20.
Byun D, Schere KL. Review of the governing equations, computational algorithms, and other components of the models-3 community multiscale air
quality (CMAQ) modeling system. Appl Mech Rev 2006;59:51–77.
Chow FK, Kosovic B, Chan S. Source inversion for contaminant plume
dispersion in urban environments using building-resolving simulations. J
Appl Meteor Climatol 2008;47:1553–72.
de Foy B, Cui YY, Schauer JJ, Janssen M, Turner JR, Wiedinmyer C. Estimating sources of elemental and organic carbon and their temporal emis25

sion patterns using a least square inverse model and hourly measurements
from the St. Louis-Midwest supersite. Atmos Chem Phys 2015;15:2405–27.
Eckhardt S, Prata AJ, Seibert P, Stebel K, Stohl A. Estimation of the vertical profile of sulfur dioxide injection into the atmosphere by a volcanic
eruption using satellite column measurements and inverse transport modeling. Atmos Chem Phys 2008;8:3881–97.
El-Harbawl M. Air quality modelling, simulation, and computational methods: a review. Environ Reviews 2013;21:149–79.
Fast J, Jr. WIG, Easter RC, Zaveri RA, Barnard JC, E. G C, Grell GA.
Evolution of ozone, particulates, and aerosol direct forcing in an urban
area using a new fully-coupled meteorology, chemistry, and aerosol model.
J Geophys Res 2006;111:D21305.
Hoek G, Krishnan RM, Beelen R, Peters A, Ostro B, Brunekreef B, Kaufman
JD. Long-term air pollution exposure and cardio-respiratory mortality: a
review. Environ Health 2013;12:43.
Hourdin F, Talagrand O, Idelkadi A. Eulerian backtracking of atmospheric
tracers. II: Numerical aspects. Q J Roy Meteor Soc 2006;132(615):585–603.
Hwang Y, Barut E, Yeo K. Statistical-physical estimation of pollution emission. Statistica Sinica 2018;in press.
Issartel JP, Sharan M, Modani M. An inversion technique to retrieve the
source of a tracer with an application to synthetic satellite measurements.
Proc R Soc A 2007;463:2863.
Karniadakis GE, Sherwin SJ. Spectral/hp Element Methods for Computational Fluid Dynamics. Oxford University Press, New York, 2005.
Keats A, Yee E, Lien FS. Bayesian inference for source determination
with applications to a complex urban environment. Atmos Environ
2007;41(3):465 –79.
Li J, Marzouk YM. Adaptive construction of surrogates for the Bayesian
solution of inverse problems. SIAM J Sci Comput 2014;36:A1163 –86.
Lieberman C, Willcox K, Ghattas O. Parameter and state model reduction for large-scale statistical inverse problems. SIAM J Sci Comput
2010;32:2523–42.

26

Martinez-Camara M, Bejar Haro B, Stohl A, Vetterli M. A robust method
for inverse transport modeling of atmospheric emissions using blind outlier
detection. Geosci Model Dev 2014;7:2303–11.
Marzouk YM, Najm HN, Rahn LA. Stochastic spectral methods for efficient
Bayesian solution of inverse problems. J Comput Phys 2007;224:560–86.
McKay M, Conover W, Beckman RJ. A comparison of three methods for
selecting values of input variables in the analysis of output from a computer
code. Technometrics 1979;21:239–45.
Phalen RF, Phalen RN. Introduction to Air Pollution Science. Jones &
Bartlett Learning, 2011.
Pope SB. Turbulent flows. Cambridge University Press, 2000.
Pudykiewicz J. Application of adjoint tracer transport equations for evaluating source parameters. Atmos Environ 1998;32:3039–50.
Rajaona H, Septier F, Armand P, Delignon Y, Olry C, Albergel A, Moussafir
J. An adaptive bayesian inference algorithm to estimate the parameters of
a hazardous atmospheric release. Atmos Environ 2015;122:748 –62.
Roosta-Khorasani F, van den Doel K, Ascher U. Stochastic algorithms for
inverse problems involving PDEs and many measurements. SIAM J Sci
Comput 2014;36:3 – 22.
Slawski M, Hein M. Sparse recovery by thresholded non-negative least
squares. In: NIPS 2011. 2011. p. 1926–34.
Stockie JM. The mathematics of atmospheric dispersion modeling. SIAM
Review 2011;53:349–72.
Thunis P, Miranda A, Baldasano J, Blond N, Douros J, Graff A, Janssen S,
Juda-Rezler K, Karvosenoja N, Maffeis G, Martilli A, Rasoloharimahefa
M, Real E, Viaene P, Volta M, White L. Overview of current regional and
local scale air quality modelling practices: Assessment and planning tools
in the EU. Environ Sci Policy 2016;65:13–21.
Tibshirani R. Regression shrinkage and selection via the lasso. J Royal
Statist Soc B 1996;58:267–88.
Tibshirani R, Saunders M, Rosset S, Zhu J, Knight K. Sparsity and smoothness via the fused lasso. J Royal Statist Soc B 2005;67(1):91–108.
27

Tibshirani. RJ, Taylor J. The solution path of the generalized lasso. Ann
Statist 2011;39(3):1335–71.
Wang M, Xu W, Tang A. A unique nonnegative solution to an underdetermined system: from vectors to matrices. IEEE Trans Signal Process
2011;59:1007–16.
Xiu D. Efficient collocational approach for parametric uncertainty analysis.
Commun Comput Phys 2007;2:293–309.
Xiu D, Hesthaven JS. High-order collocation methods for differential equations with random inputs. SIAM J Sci Comput 2005;27:1118–39.
Xiu D, Karniadakis G. The Wiener-Askey polynomial chaos for stochastic
differential equations. SIAM J Sci Comput 2002;24:619–44.

28

