SCARLET: Source separation in multi-band images by Constrained Matrix Factorization
Peter Melchior Fred Moolekamp Maximilian Jerdee Robert Armstrong Ai-Lei Sun James Bosch Robert Lupton




1]Peter Melchior
peter.melchior@princeton.edu
1]Fred Moolekamp
fredem@princeton.edu
1]Maximilian Jerdee
1]Robert Armstrong
[1]Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA
2]Ai-Lei Sun
[2]Department of Physics & Astronomy, Johns Hopkins University, Bloomberg Center, 3400 N. Charles St., Baltimore, MD 21218, USA
1]James Bosch
1]Robert Lupton


We present the source separation framework  for multi-band images, which is based on a generalization of the Non-negative Matrix Factorization to alternative and several simultaneous constraints.
Our approach describes the observed scene as a mixture of components with compact spatial support and uniform spectra over their support. 
We present the algorithm to perform the matrix factorization and introduce constraints that are useful for optical images of stars and distinct stellar populations in galaxies, in particular symmetry and monotonicity with respect to the source peak position.
We also derive the treatment of correlated noise and convolutions with band-dependent point spread functions, rendering our approach applicable to coadded images observed under variable seeing conditions.
 thus yields a PSF-matched photometry measurement with an optimally chosen weight function given by the mean morphology in all available bands.
We demonstrate the performance of  for deblending crowded extragalactic scenes and on an AGN jet – host  galaxy separation problem in deep 5-band imaging from the Hyper Suprime-Cam Stategic Survey Program.
Using simulations with prominent crowding we show that  yields superior results to the HSC-SDSS deblender for the recovery of total fluxes, colors, and morphologies.
Due to its non-parametric nature, a conceptual limitation of  is its sensitivity to undetected sources or multiple stellar population within detected sources, but an iterative strategy that adds components at the location of significant residuals appears promising.
The code is implemented in Python with C++ extensions and is available at <https://github.com/fred3m/scarlet>.






 methods: data analysis techniques: image processing galaxies: structure Non-negative matrix factorization









§ INTRODUCTION



Modern astronomical wide-field surveys cover large areas of the sky at ever increasing depths, revealing more objects and low-surface-brightness features of extended objects that were previously too faint to detect.
These gains drive investigations into galactic, extragalactic and cosmological phenomena at an unprecedented level of detail and statistical power.
On the other hand, because of the enhanced sensitivity, a larger fraction of the observed area is associated with detectable objects, thereby increasing the chance that multiple objects overlap.
This so-called “blending” constitutes a major concern for the analysis of existing and upcoming surveys, especially those that observe from the ground.

The majority of methods for measuring the properties of celestial objects assume that every object can be considered isolated.
If that assumption holds, well-defined and accurate measurements of the flux, position, shape, and morphology can routinely be made with methods that are either based on moments of the light distribution within some aperture or on parametric fits to the images. 

However, the notion of isolated objects is becoming increasingly obsolete.
With a limiting magnitude of i≈24, DES[<https://www.darkenergysurvey.org>] <cit.> finds that 30% of galaxies suitable for weak-lensing measurements are affected by blending <cit.>. 
For HSC[<http://hsc.mtk.nao.ac.jp/ssp/>] <cit.>, whose Wide survey has a limiting magnitude of i≈26, <cit.> find that 58% of measured objects are in blended groups, a dramatic increase despite a substantially better average seeing than DES.
LSST[<http://lsst.org>] <cit.> expects to reach i≈27 after 10 years of operations , and it is estimated that 63% of observed galaxies will have Sérsic model photometry that is altered by more than 2% due to the presence neighbors  (Sanchez et al., in prep.).
Even more problematically, in a comparison study of HST and Subaru imaging of a galaxy cluster field, where the Subaru data had similar seeing and depth so as to serve as a proxy for LSST, <cit.> found that 14% of observed galaxies are blended but not recognized as such in the ground-based images.
Slightly larger numbers for unrecognized blends are found for HSC in a study that inserted fake objects into  real survey images to infer how many of them could be recovered (Murata et al., in prep.).

This lack of separability between objects necessitates the employment of techniques that analyze entire scenes with overlapping objects.
For direct measurements, such as moments within apertures, there is no accurate way to correct for the excess light from overlapping objects because such a correction depends e.g. on radial profiles of the objects involved, which cannot be determined well for blends.
As a consequence, multiple objects need to be modeled either iteratively (by masking all but one) or simultaneously, often requiring sophisticated and fine-tuned schemes to prevent unstable or physically implausible solutions <cit.>.
Any such scheme is suitable to extract and separate the objects in celestial scenes, i.e. to “deblend” those scenes.
Differences between the schemes include the propagation of errors, which conceptually favors simultaneous approaches, and whether the desired measurements are generated directly from deblender models or by passing them on to established measurement algorithms for isolated objects.

Traditional deblending approaches in astronomy are achromatic, i.e. they employ information from only a single image.
SExtractor detects blending by thresholding an image at a range of intensity levels and searching for sets of pixels that are connected at a lower threshold but split into several connected regions at a higher one <cit.>.
As a consequence of the splitting approach, the association of pixels to objects is unique and exclusive, i.e. in the internal representation of blended objects they do not overlap.
This unrealistic notion has necessitated mitigation strategies or fine-tuning to prevent “over-deblending” of larger galaxies caused by smaller and fainter companions or interlopers <cit.>.

The deblender in the SDSS Photo pipeline (Lupton, in prep.) does allow for overlap between nearby objects and, consequently, needs to estimate the portion of any pixel's flux that is due to each object.
It uses a two-step approach, in which first a template is constructed for each object based on the requirement that pixel values symmetrically across the object's peak pixel be identical; they are generally not, so the minimum pixel value of those two pixels is adopted for both.
Then, the original image values are projected onto those templates, associating each pixel's flux to different objects in proportion to the amplitude of the respective templates.




Despite very few assumptions, the method mostly separates sources into physically plausible objects but struggles with situations where a central object is symmetrically surrounded by neighbors, for instance a blend with three peaks in a row (Lupton, in prep.).



More recently, <cit.> and <cit.> proposed variants of inpainting techniques, where the relevant pixels of blended objects are replaced by an estimate of a local variable background, working inwards from an initially defined outline.
The portion of the pixel flux above the background estimate is attributed to the respective object.
These approaches implicitly account for blending by assuming that the background captures the flux contributions from neighboring objects.
This is particularly useful when recovering small objects in multi-scale blending situations like the cores of galaxy clusters and  removing them from the scene so that large objects can be measured separately.

While effective in many cases, all of the deblending schemes outlined above employ heuristic arguments for how to separate overlapping sources.
They also perform the pixel–object association sequentially, one object at a time, thus losing the advantages of a simultaneous solution, for instance the ability to explore the degeneracies that arise because the objects are not isolated. 
However, it is our opinion that the biggest limitation stems from the restriction to a single image, and therefore a single filter band, while most modern surveys observe the sky in several filters. 
A visual inspection of multi-band images clearly suggests that color can serve as a powerful discriminator between different objects, even with severe overlap (see <ref>).

MuSCADeT <cit.> addresses both limitations by building a joint model of multi-band image data.
As their model is non-parametric, the number of degrees of freedom is large, which leads to many possible degeneracies in the solutions, so they demand that the spatial distribution of each source be sparse in the starlet (a form of isotropic wavelet) domain.  
The resulting solutions extract preferentially compact features down to the noise level, using a set of previously identified colors for each feature.
For applications to wide-field multi-band data, we cannot generally assume to a priori identify the color of an object in the scene because there might not be a single pixel whose color is uncontaminated by other objects.
We therefore seek the ability to update both spectral and morphological characteristics of the objects.
While sharing noticeable similarities with MuSCADeT regarding the use of a non-parametric constrained morphological model, one can consider our approach an extension that also updates the source spectra as well as a generalization that allows an arbitrary number of constraints to be placed on each source.

The outline of the paper is as follows: We introduce our approach, dubbed , in <ref>, demonstrate its performance on real data and simulations in <ref>, and conclude in <ref>. 



§ METHODOLOGY



We base our deblender on the assumption that astronomical scenes are superpositions of multiple components, each with


  * a spatially compact support and

  * a constant spectrum over that support.

For stellar fields this is obviously true, but even for galaxies, especially marginally resolved ones, which constitute the vast majority of galaxies in deep surveys, the assumption is appropriate.
In addition, even large, extended galaxies can be thought of as conglomerations of components (e.g. bulges, discs, bars, star forming regions) for which the assumptions above hold at least approximately.
[In the literature the terms “source separation” and “component separation” are often used interchangeably. To better reflect the hierarchical nature of astrophysical scenes, we will define the term “Source” as a collection of co-centered components that belong to the same astrophysical object.]
For instance, the popular bulge-disc decomposition for galaxies <cit.> is justified by this interpretation.
We note that the assumption of linear superposition implies that components do not interact, which is correct only for transparent emitters.
Absorption, e.g. by dust in the galaxy, can be approximated by allowing negative values in the source spectrum, but substantial opacities cannot fully be modeled because the effect depends on the amount of absorbing material and the intensity of the background radiation.

The assumptions above appear to lend themselves to a parametric modeling framework, where one assumes to know the shapes of the components and potentially their intrinsic spectrum, exploiting quite tight relations between colors and morphologies exhibited by galaxies in the late universe <cit.>.
While drastically reducing the number of optimization parameters, we are critical of this approach for two reasons:
First, in the translation of an intrinsic, restframe spectrum to the observed broadband colors one needs to take the galaxy's redshift into account, which is equivalent to estimating a photometric redshift as part of the deblending process.
If the redshifting prescription is incorrect, e.g. because of a limited library of spectra or the evolution of those spectra with redshift,
it would affect the properties of the deblended components—not only their recovered spectra but also their shapes.
Second, at the stage in the analysis pipelines of large astronomical surveys where we envision the deblender to operate, it will not necessarily be established what kind of sources are in the scene; in other words, a suitable parameterization is probably not known.
This is most evident when looking at the star-galaxy distinction: for stars, three parameters are sufficient (two centroid coordinates and one amplitude), while even simple galaxy models need at least one more parameter (the size).
Model-fitting under those conditions can be done by transdimensional sampling <cit.>, but the computational costs are likely too high for large-volume data sets.

Because of these concerns, we seek to characterize the scenes without making questionable astrophysical assumptions, which means describing colors in the observed frame and morphologies in the free-form space of image pixels.



 §.§ Non-negative Matrix Factorization



We assume that an astronomical scene  that we seek to analyze is organized in the form of a multi-band image cube of aligned images in B bands, each of which is suitably flattened to have a total number of N pixels.
Our previous assumptions give rise to a multi-band model  as a sum of a finite number of components K,

    = ∑_k=1^K ^⊤_k ×_k = ,

where _k∈ℝ^B is the amplitude of component k across all bands, i.e. its spectral energy distribution (SED), and _k∈ℝ^N is the spatial shape of that component. By arranging the _k as columns of ∈ℝ^B× K and _k as the rows of ∈ℝ^K× N, we have a model  that is given by the product of two matrix factors.

With a homoscedastic Gaussian error model, which is appropriate for most extragalactic images in the optical, the likelihood function is

    f(,) = 1/2-_2^2

where ._2 denotes the element-wise L_2 (Frobenius) norm.
[The objective function in <ref> is insufficient for dealing with realistic multi-band data, which generally exhibit correlated noise and correlated signals, the former from warping the images to rectify any astrometric distortion, the latter because of the blurring from the point spread function.
We will work out in <ref> and <ref> how to deal with both effects.]
In its simplest form, the Non-negative Matrix Factorization <cit.> then amounts to fitting  and  such that they minimize f and obey the non-negativity constraint, which is given by the indicator function of the  set of non-negative matrix elements:

    g_+(X) = 
    0    if X_mn≥ 0 ∀ m,n
    ∞   else.

In other words, one seeks to minimize f(,) + g_+() + g_+().
The classical way of solving this NMF problem is known as “multiplicative updates” <cit.>, which suffers from poor convergence if the constraints strongly work against the minimum of the objective function.
Moreover, in its simplest form the NMF is hampered by a degeneracy that stems from the transformation (,)→(Q,Q^-1) with an arbitrary invertible matrix Q, which means power can be shifted from  to  or vice versa, hampering convergence to a well-defined solution.
To compensate one can introduce another constraint on the norms of either  or .
The most natural normalization for astronomical scenes is to require that each component's SED sums up to unity.
We therefore seek to minimize f(,) + g_+() + g_unity() + g_+(), where g_unity is another indicator function of the set of matrices in which the columns are properly normalized.

We do so by using a two-block variant of the “proximal gradient method”.
In short, a closed smooth proper convex function f(x⃗) under a convex constraint function g(x⃗) with x⃗∈ℝ^n can be minimized with the following sequence:

    x⃗^it+1←_λ^it g(x⃗^it - λ^it∇ f(x⃗^it)).

The step size λ must be chosen from (0,2/L], where L is the Lipschitz constant of ∇ f.
The function g is not directly evaluated but accessed through its proximal operator,

    _λ g(x⃗)≡v⃗argmin{ g(v⃗)+1/2λx⃗-v⃗_2^2}.

The proximal operator of g effectively hides the infinite values and the non-differentiability of its underlying function.
For instance, the proximal operator of the indicator function I_𝒞 of a non-empty set 𝒞 is the Euclidean projection operator onto the set—i.e. it yields v⃗∈𝒞 that is nearest to —and the scaling parameter λ is irrelevant.
In particular, 

    _+(x⃗) = (x_1^+, …, x_n^+) with x_i^+ ≡max(0,x_i)
       _unity(x⃗) = x⃗_1^-1x⃗,

For many other classes of constraint functions proximal operators have analytic forms, which means that the minimization of <ref> does not have to be carried out.
For more details on proximal methods, we refer the reader to <cit.> and <cit.>.

As long as there is only one constraint per variable, we can interleave update steps of the form of <ref> that minimize f( | ) and those that minimize f( | ), effectively a proximal implementation of a coordinate descent method, as shown in <ref>. Note that we merged g_+() + g_unity() into one proximal operator _unity+≡_+(_unity), which can be done because the two constraints commute, i.e. they do not interfere.
This reduces the problem to calculating the gradients

    ∇_ f(,) = -(-) ^⊤
       ∇_ f(,) = -^⊤(-)

and their Lipschitz constants

    L_ = ^⊤_s and L_ = ^⊤_s,

where ._s refers to the spectral norm. 
The sequence in <ref> approaches a fixed point that is also a local minimum of f and can therefore be terminated whenever the change between successive iterations falls below a user-defined threshold.

This approach is similar to the one of <cit.> in that the minimization is performed in two blocks and that constraints are expressed by their proximal operators, but we only perform one gradient step instead of seeking the exact solution for every block in every iteration.
This reduces the number of gradient and constraint evaluations while still yielding a sequence that is guaranteed to converge <cit.>.

 



 §.§ Constrained Matrix Factorization


For many practical applications we want to be able to impose additional constraints besides (or instead of)  non-negativity.
We will give a list of alternative constraints in <ref>, but it is useful to think of constraints more generally as logarithms of priors on the solution.
In addition, we are often forced to impose multiple simultaneous constraints to deal with the large number of unknowns in <ref>, especially when the number of components K is larger than the number of bands B.
To reflect these more general conditions, we favor the term “Constrained Matrix Factorization” (CMF) instead of NMF.

With the direct projection method of <ref>, we can only accommodate any one proximal operator per variable, even though non-interacting constraints can be daisy-chained to act like a single one, as we have done with _unity+.
One could generalize this approach by the method of Alternating Projections, also known as Projection onto Convex Sets <cit.>, but that still requires that the proximal operators act directly on the variables  and .
Often, it is beneficial to express constraints in a transformed domain, for instance demand that the solution be sparse in that domain <cit.>.
We therefore prefer the Alternating Direction Method of Multipliers <cit.>.
In its basic form it seeks to

    x⃗minimize f(x⃗)+g(Lx⃗)

which allows us to introduce an arbitrary matrix L.
For this work L will be a gradient or symmetry operator, but it could as well be a Fourier or wavelet transform, etc.
The previous equation can be re-written in “consensus form”

    [    minimize f(x⃗)+g(z⃗);  subject to   Lx⃗-z⃗=0, ]

which suggests an approach that splits the optimization into two separate tasks: one that minimizes f and another that satisfies g by introducing the auxiliary variable z⃗ <cit.>.
The updating scheme also introduces the dual variable u⃗, which connects the two tasks, creating the following sequence:

    [                         ^it+1                             ← _λ f(^it-λ/ρ^⊤(^it-^it+^̆it));                         ^it+1                             ←              _ρ g(^it+1+^̆it);                        ^̆it+1                             ←             ^̆it+^it+1-^it+1, ]

where the scaling parameters need to satisfy 0<λ≤ρ/||||_s^2.
[For the -update, we use a linearized form of the ADMM, which is known as split inexact Uzawa method <cit.>.]
Two aspects of this sequence are remarkable. 
The function f is only accessed through its proximal operator, which means it does not need to be differentiable.
In our case, f is differentiable, and we can consider

    _λ f() ≡x⃗ - λ∇ f()

the corresponding proximal operator of the first-order approximation to f.
[In fact, we can treat the entire projected gradient step of <ref> as a first-order approximation of f + g, which means we can satisfy a single constraint in the direct domain, such as non-negativity, without having to invoke the auxiliary variables  and $̆.
Second, the linear operatordoes not need to be invertible; it and its transpose are only used for dot products, rendering the updates particularly efficient with sparse matrices. 

In <cit.>, we have demonstrated that the ADMM can be extended to multiple variables of a functionfthat is convex in each of those variables by constructing a block-updating sequence analogous to <ref>.
For the CMF,and(which we will call_1and_2and enumerate with the indexjfor the remainder of this subsection) thus have their own set of ADMM auxiliary variables_1, _1,_2, _2. 
In addition, we have introduced a method that allows each_jto be subject to multiple simultaneous constraints, one pair (_ji, _ji)for each constraintg_ji(_j)with i∈{1,…,M_j }.
Combined, our new method, dubbed , provides the flexibility we seek for deblending astronomical scenes via the CMF.
It will 

    [                           _1, _2minimize   f(_1, _2) + ∑_j=1^2 ∑_i=1^M_jg_ji(_ji);                               subject to _ji_j-_ji=0  ∀ j∈{1,2} and i∈{ 1,…,M_j } ]

through the Block-Simultaneous Method of Multipliers <cit.>.

The core working principle of  is given in <ref>.
The ADMM updates of <ref> are performed in lines 10–13, using the gradient updates from <ref> in the form of <ref> in line 10.
The step sizesλ_jin line 8 are given as before by1/L_jwith the Lipschitz constants from <ref>.
It is worthwhile pointing out that in the limit of all constraints being satisfied,=and=0, in which case <ref> is equivalent to <ref>.

Following <cit.>, we implement stopping criteria (on line 14) based on primal residualP^it+1=^it+1-^it+1and the dual residualD^it+1=1/ρ^⊤(^it+1-^it).
These residuals describe how close the currentis to satisfying the constraint in the transformed domain and how muchchanges, respectively.
To assess primal and dual feasibility, we require

    P^it+1_2 ≤ϵ^pri≡√(p) ϵ^abs + ϵ^relmax{L^it+1_2, ^it+1_2}  and
       D^it+1_2 ≤ϵ^dual≡√(n) ϵ^abs + ϵ^rel / ρ L^⊤^it+1_2,

wherepandnare the number of elements inand, respectively.
The error thresholdsϵ^absandϵ^relcan be set at suitable values, depending on the precision and runtime constraints of the application.

The procedure outlined in <ref> yields point estimates of the matrix factors.
Uncertainties on each factor are currently calculated by linear error propagation based on <ref>.
While those uncertainties provide estimates of the significance of each element of the matrix factors, they do not properly capture the interdependence of the model components or the effects of the constraints on the degeneracies between them.
Fully accurate uncertainties are difficult to estimate because of the non-differentiability of the constraint functions.
We plan to investigate ways to better characterize the uncertainties of the  model in forthcoming work.

This concludes our presentation of the general approach we have adopted in .
A reader who is mostly interested in the astrophysical applications can safely jump to <ref> and skip details of the formalism below.

[t]
Block-SDMM Constrained Matrix Factorization The matrix factors  and  are treated as two blocks X_j with j∈{1,2} and independently updated with a sequence of gradient steps _λ_j f,j and jointly optimized with an arbitrary number of constraint functions g_ji with i∈{1,…,M_j}.
All constraints are expressed as proximal operators _ρ_ji g_ji and include transformation matrices _ji (cf. <ref> for the generic ADMM update sequence). The parameters λ_j and ρ_ji are the maximum allowable step sizes in this scheme.[1]
CMF-BSDMM^0, ^0_1^0, ^0_2 ^0, ^0j=1,2_ji^0 _ji_j ∀ i ∈{ 1,…,M_j}_ji^0 0⃗ ∀ i ∈{ 1,…,M_j}it = 0,1,…j=1,2λ_j^it j==1? _2^it_2^it⊤_s^-1 :  _1^it+1⊤_1^it+1_s^-1ρ_ji^it 2 M_j λ_j^it+1_ji_2^2_j^it+1_λ_j f,j(_j^it-∑_i=1^M_jλ_j/ρ_ji_ji^⊤(_ji_j^it-_ji^it+_ji^it))i = 1,…,M_j_ji^it+1_ρ_ji g_ji(_ji_j^it+1+_ji^it)_ji^it+1_ji^it+_ji_j^it+1-_ji^it+1⋀_ji{P_ji^it+1_2 ≤ϵ^pri∧D_ji^it+1_2 ≤ϵ^dual}
     break
    

 §.§ Constraints


The conventional NMF is capable of separating components, but only from the color contrast between objects in the scene.
If the color contrast vanishes or the numberKof overlapping components gets large, the model in <ref> becomes increasingly susceptible to degeneracies, where the mixture model yields a fair representation of the scene but many of the individual components look implausible.

We thus prefer the flexibility of the CMF to impose additional priors or constraints based on existing knowledge or intuition about how SEDs and morphologies are distributed.
In addition to non-negativity, which is justified for astronomical sources that emit photons, we will routinely employ symmetry, monotonicity, and sparsity constraints for the galaxy morphologies (see <ref>).
We currently do not utilize priors on galaxy or stellar SEDs, but we will demonstrate in <ref> how color constraints can help in cases when additional stability of the source extraction is desired.

We want to emphasize that  is capable of incorporating any prior or constraint that can be expressed as a proximal operator.
We therefore think of it as a general-purpose component separation framework that can be optimized for a wide variety of applications.



  §.§.§ Symmetry

We take guidance from the SDSS  Photo deblender (Lupton, in prep.), which, inspired by <cit.>, constructs templates of the sources that are symmetric under180^∘rotations around the peak pixel.
We formalize this requirement by stating that 

    ∀ k,i: _k,i - _k,j = 0,

whereiandjare the indices of the symmetric partner pixels,
[In order to always have a well-defined partner for each pixel (other than the peak pixel), the peak pixel needs to be in the center of the image, and the number of pixels in the image must be odd in both directions.
We achieve this by defining a virtual box around each source that contains all of its pixels, with the peak pixel at the box center.

]
which can be expressed in matrix form,

    ∀ k: _symm_k = 0⃗,

with a symmetry matrix_symm∈ℝ^(N/2-1) ×Nthat has the elements(1,-1)for each ofN/2-1suitable pixel pairs (the peak pixel is excluded).
In this linearly transformed space, the constraint is imposed by the projection operator onto the zero vector, i.e._0(_symm_k).
The combination of the transformation matrix and the proximal operator is required for any ADMM-style constraint in <ref>.

Compared to the SDSS  Photo deblender, which adopts the minimum value for both pixels in the symmetric pair, this operator will converge to the mean of both pixels.
This is a direct consequence of the proximal operator being the Euclidean projection operator and the mean statistic having the least squared error.



  §.§.§ Monotonicity

The prominent failure mode of the SDSS  Photo deblender arises from its templates not necessarily declining monotonically from the peak pixel.
[To address this weakness, the implementation of the SDSS  Photo deblender in the current HSC software stack imposes a monotonicity constraint on previously smoothed templates.]
While resolved images of real galaxies may in detail show non-monotonic behavior, for instance due to star-forming regions in the outer parts of spiral galaxies, we prefer imposing this constraint for the bulk source shape.
An area that violates this constraint likely constitutes a different stellar population, potentially with a color different from the bulk model, and should thus be fit by adding another (monotonic) component to the mixture model.

Monotonicity for two-dimensional features is not uniquely defined.
We use a form of radial monotonicity, i.e.

    ∀ k,i: ∑_j∈𝒰(i) w_ji_k,j - _k,i≥ 0,
  
where𝒰(i)is the set of neighbor pixels ofithat are closer to the peak thani, and the weights are normalized such that∀i: ∑_j w_ji=1.
We implement two variants:

  *  The neighborhood contains only the nearest neighbor of i in the direction of the peak, 𝒰(i) = {NN(i) } and w_NN(i),i = 1.

  *  The neighborhood contains all three neighbor pixels that are closer to the peak, and their weights are given by w_ji = cos(ϕ_j - ϕ_i), where ϕ_i denotes the angle from the center of pixel i to the center of the peak pixel.

Unsurprisingly, we find the latter variant to yield models that are somewhat softened in azimuthal direction, and they can violate monotonicity along the straight path from pixelito the peak, which variant 1 would obey.

As before, we rewrite the constraint in matrix form,∀k: _mono _k ≥0, which can be enforced by the non-negativity constraint from <ref> in the transformed domain,_+(_mono ).

This enforcement of monotonicity can be slow to converge in heavily blended regions.
The   package therefore contains a projection operator in the direct domain,_mono,  which constructs a monotonic version of its argument. While this construction will in general not yield the  monotonic solution that is closest to the argument (the operator is thus in fact not a true proximal operator), we have found it to be substantially more robust and efficient in practice.



  §.§.§ Sparsity


Sparsity refers to a solution that can be represented with a small number or a low amplitude of coefficients in a particular domain.
In the direct, i.e. the pixel, domain, sparsity constraints thus induce solutions with few non-zero pixels and/or compressed pixel amplitudes, which can be useful to reduce the impact of pixel noise.
They are introduced by minimizing the objective function and a penalty function based on the norm of the solution, often theℓ_2orℓ_1norms or theℓ_0pseudo-norm, or a combination of them.
While theℓ_2norm is differentiable, the other two are not. However, bothℓ_0andℓ_1have simple proximal operators, known as “hard thresholding” and “soft thresholding” operators, respectively:

    _λℓ_0() = (x_1^h, …, x_n^h) with x_i^h = 0    if |x_i|<λ
     x_i    else
       _λℓ_1() =  (x_1^s, …, x_n^s) with x_i^s = sign(x_i) max(0, |x_i|-λ)

This is the first example of a proximal operator of a function other than an indicator of a closed set.
The parameterλthen governs the strength of the penalty compared to the objective functionf.

Because  allows arbitrary transformations ofandwith linear operators, the sparsity constraints can be applied in the direct domain and/or in the transformed domain. 



  §.§.§ Other options


We do not attempt to provide a comprehensive list of possible constraints that are expressible by proximal operators. Instead, we discuss several that we found useful and refer the reader to the appendix of <cit.> and section 6 of <cit.> for more details.

Spatially flat components, e.g. for fitting a uniform sky emission, can be realized by projecting onto the mean (see discussion in <ref>):_flat() = 1n∑_i^n x_i 1⃗_nwith∈ℝ^n. 

To allow some amount of spatial variability, e.g. for describing Intra-Cluster Light, we find the Maximum Entropy regularization <cit.> useful, which fundamentally seeks to find the flattest distribution that is compatible with the data by adding the (negative) Shannon entropy term

    S = ∑_i^N p_i log(p_i) with p_i = x_i / ∑ x_i

to the objective functionf.
The proximal operator ofSis 

    prox_λ S() = λ[W(1/λexp(/λ-1))],

whereWdenotes the Lambert-Wfunction.






By choosing the matrixto represent a basis transform, one can fundamentally alter the characteristics of the solution one wants to isolate or penalize.
A few examples with applications in data analysis are the Fourier transform <cit.>, the wavelet transform <cit.>, or the shapelet transform <cit.>.
In the transformed domain, simple projection operators can then e.g. null undesired frequencies.

Finally, instead of hard constraints one can also introduce prior distributions on the solutions.
If the prior distributionp()is smooth, the corresponding proximal operator is given by the gradient step oflog(p)analogous to <ref>; if it has hard edges it can be represented by wrapping the gradient step in a proximal gradient update as in <ref>.



 §.§ Heteroscedastic and correlated errors
<ref> does not account for heteroscedastic errors.
<cit.> introduced a simple extensions to heteroscedastic errors with the matrixVof variances of each pixel in:

    f(A,S)=V^-1/2∘(A·S - D)_2^2,

where∘denotes the element-wise, or Hadamard, product.
It is straightforward to adjust the gradients in <ref> but it is not clear how their Lipschitz constants differ from <ref>.
[The two works cited above did not need to compute the respective Lipschitz constants as they solved the NMF by multiplicative updates, not proximal algorithms.]
In addition, <ref> does not account for correlations in the noise that arise when images are astrometrically rectified.

For that purpose it is necessary to serialize,, andinto vectors (Y⃗,A⃗,S⃗) with lengthsB·N,B·KandK·N, respectively.
The most convenient vector formY⃗ = (_1, …, _B)simply concatenates the images inBbands.
The covariance matrix ofY⃗has block-diagonal form

    Σ_Y⃗ = [ Σ_1   0   …   0;   0 Σ_2   …   0;   ⋮   ⋮   ⋱   0;   0   0   … Σ_B;     ],

withΣ_b∈ℝ^N×Ndenoting the pixel covariance matrix in bandb.
We can now rewrite the objective function in a form that uses vectors forand one ofor, with the other one expressed by a sparse block matrix (or).
The objective function is then

    f(A⃗,) = 1/2(Y⃗-A⃗)^⊤Σ^-1_Y⃗ (Y⃗-A⃗)_2^2,
       f(,S⃗) = 1/2(Y⃗-S⃗)^⊤Σ^-1_Y⃗ (Y⃗-S⃗)_2^2

which describes an ordinary vector-valued linear inverse problem.
The gradients are

    ∇_A⃗ f(A⃗,)    = -^⊤Σ_Y⃗^-1(Y⃗-A⃗)
    ∇_S⃗ f(,S⃗)    = -^⊤Σ_Y⃗^-1(Y⃗-S⃗)

and the covariance matrices are given by

    Σ_A⃗   = ^⊤Σ_Y⃗^-1
    Σ_S⃗   = ^⊤Σ_Y⃗^-1.

The Lipschitz constants of those gradients, which we need to calculateλ_jin <ref>, are the spectral norms of the covariance matrices.
Finally, the block matrices that encodeorare given by

    = [ _11I_N _21I_N      … _K1I_N; _12I_N _22I_N      … _K2I_N;      ⋮      ⋮      ⋱      ⋮; _1BI_N _2BI_N      … _KBI_N;        ]∈ℝ^BN× KN

and 

    = [ ^⊤  0  …  0;  0 ^⊤  …  0;  ⋮  ⋮  ⋱  0;  0  0  … ^⊤;    ]∈ℝ^BN× BK,

each stacking the blocks row-wise over multiple bands.



 §.§ PSF (de-)convolution


Incorporating the convolution with the point spread function (PSF) into our model is particularly useful for a deblender because it allows a distinction between overlap that is caused by the physical extent of the sources involved and the one due to atmospheric or instrumental blurring of the images.
Knowing the latter reduces the ambiguity caused by the former.

In addition, our model from <ref> fits the same morphology to each band, which
is only meaningful if the PSF does not change between bands or the images have been PSF-homogenized across all bands.
Otherwise, the data will exhibit color variations that our model cannot reproduce, even if the underlying morphological structure of the astronomical source is achromatic like e.g. that of a star.

In the vectorized frame we have introduced in <ref> we express the convolution as

    P̅ = [ P_1   0   …   0;   0 P_2   …   0;   ⋮   ⋮   ⋱   0;   0   0   … P_B;     ]∈ℝ^BN× BN,

with each of the blocks describing the PSF in a given band as a linear operator that couples pixels ofS⃗.
The model can thus be extended to have the objective functions

    f(A⃗,) = 1/2(Y⃗-P̅A⃗)^⊤Σ^-1_Y⃗ (Y⃗-P̅A⃗)_2^2,
       f(,S⃗) = 1/2(Y⃗-P̅S⃗)^⊤Σ^-1_Y⃗ (Y⃗-P̅S⃗)_2^2

with gradients

    ∇_A⃗ f(A⃗,)    = -(P̅)^⊤Σ_Y⃗^-1(Y⃗-P̅A⃗)
    ∇_S⃗ f(,S⃗)    = -(P̅)^⊤Σ_Y⃗^-1(Y⃗-P̅S⃗)

and covariance matrices

    Σ_A⃗   = (P̅)^⊤Σ_Y⃗^-1P̅
    Σ_S⃗   = (P̅)^⊤Σ_Y⃗^-1P̅,

from which we can compute_fand the step sizesλ.

This ansatz with the PSF as a linear operator is only valid as long as the PSF images and the deconvolved model remain well sampled.
[Internal oversampling is an option to relax this requirement, which we will investigate in future work.]
Otherwise, the discretized nature of either will lead to convolution artifacts, which may result, e.g., in apparent non-monotonic behavior in the deconvolved morphology even if the source is perfectly monotonic.

We therefore recommend to construct difference kernels for each band,P_Δ,b P_min = P_b, such as to perform a partial deconvolution, resulting in an image with a band-independent PSFP_min.
Requiring Nyquist-sampling forP_Δ,bandP_minsuggest thatP_minshould have a widthTof about12 min_b T(P_b).



 §.§ Translations and centering


In blends, peak positions from the initial detection step may not accurately reflect the true source centers.
In fact, the apparent peak positions are often shifted towards neighboring objects.
In addition, the combination of symmetry and monotonicity constraints results in the model being strongly affected by the assumed location of each source's center.
To avoid spurious dipole residuals from inaccurate center positions, we introduce linear translation operatorsT_xandT_y∈ℝ^N×Nto shift the position of each source.
These operators encode the translation operation, e.g.

    ∀ i<N: (_k,i, S_k,i+1) → ((1-dx) S_k,i, dx S_k,i+1)

for a shift of componentkby an amountdxin the positive horizontal direction, as sparse band-diagonal matrices, i.e. we treat the translation as a linear interpolation between the pixel and its relevant neighbors.
This is again feasible only if the source is well sampled, and higher-order interpolation schemes can be implemented if needed.

We introduce the translation in horizontal and vertical directions into the model via𝐒^'= T̅ 𝐒, using the block-diagonal matrix

    T̅ = [ T_1   0   …   0;   0 T_2   …   0;   ⋮   ⋮   ⋱   0;   0   0   … T_K;     ]∈ℝ^KN× KN,

where the two-dimensional shift for each component is given byT_k = T_y,k T_x,k.
If multiple components belong to the same source, they have identicalT_k.
The shifted morphology vector𝐒^'is then used instead of𝐒in <ref>, <ref>, and <ref>.
This way, all gradients and Lipschitz constants are correctly computed with the equations given above.

To determine optimal source centers, we create a new set ofT^'_x,kandT^'_y,kmatrices, each with an additional shiftδof typically0.1 pix in either horizontal or vertical direction.
We form two difference images for each component,

    D_x,k   = A_k  ×( T^'_x,k - T_k) S_k
    D_y,k   = A_k  ×( T^'_y,k - T_k) S_k,

and combine them into the matrix

    D̅ =
        [ D_x,1 D_y,1 D_x,2 D_y,2     … D_x,K D_y,K ].

We then state that the current model residuals are a linear combination of the difference images,

    Y⃗-P̅S⃗ = D̅Δ,

an ordinary linear inverse problem for the two-dimensional offsetsΔ, which we solve with a least-squares solver.
We perform this secondary minimization every 10-th iteration, so that the updates forhave time to adjust to changes in the centers.



§ APPLICATIONS


 §.§ Scene from the HSC UltraDeep survey

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >


    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >


    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >


    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >


    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

Deblender comparison for several sources (one per row) from <ref>. 1st column: HSC grizy coadd data. 2nd column: The corresponding component of the  model. 3rd column: Data from the 1st column, restricted to the object's pixels in the  SExtractor segmentation map of the detection image. 4th column: HSC-SDSS deblender result for the object, run independently on each band.

We demonstrate the performance of  in deblending deep ground-based images from HSC.  
The test case shown in <ref> stems from the UltraDeep COSMOS coadd images of the HSC public data release <cit.>.
The data are taken in filters  grizy to a 5σpoint source depth of 27.4, 27.3, 27.0, 26.4, and 25.6, respectively, with average seeing of 0.74 arcsec or better.
This example was selected because of its complex configuration with 33 detected sources, which span a wide range of morphologies, sizes, and brightnesses.

We detect the sources in a detection image, defined as the sum across all bands, using the  Python package  sep[<https://github.com/kbarbary/sep>]<cit.>, which implements the functionality of  SExtractor.
The detection threshold was set at1.2σabove the background, with very sensitive deblending parameters  DEBLEND_THRESH=64 and  DEBLEND_CONT=3e-4; these settings were chosen to detect the majority of sources in this particular example, without over-deblending the larger galaxies in the image.



  §.§.§ Results of 

We specify the CMF model in the simplest form, namely with single-component sources, which means the model cannot describe color variations within a galaxy.
We initialize the component morphologies_kby constructing symmetric, monotonic approximations to the region in the detection image around each identified peak, very similar to the approach the SDSS  Photo deblender employs.
The initial SEDs_kare set to the color of the data averaged with the spatial weight function_k.
We perform 200 iterations of <ref>, using inverse per-pixel variances as weights, but ignoring any noise correlation from resampling on the coadd image.
We also ignore the relatively small PSF variation between bands in the example, thereby performing the source separation in the convolved frame.

As one can see in the central panel of <ref>, most features of the data are represented well by the  model, despite the limited flexibility of one-component sources.
Furthermore, there are two resolved spiral galaxies (IDs 19 and 22), which in detail violate our assumptions of symmetric and monotonic behavior.
Nonetheless, the residuals in the right panel of <ref> indicate that this simple model can describe most of the source morphologies.
Moreover, the residuals clearly reveal that our model could be extended by adding several components to account for previously undetected localized emission, either from unrelated sources (such as the red clumps in the regions of IDs 19 and 25) or distinct stellar populations in detected sources (as in the centers of IDs 24, 29, and 33).

We leave the model refinement based on the analysis of the residuals to future work, and instead compare the deblending results of , the HSC-SDSS deblender, and  SExtractor for several sources in <ref>.
Because  SExtractor does not attempt to separate overlapping sources (it merely seeks to find a suitable boundary between them), we show the multi-band data for the pixels belonging to the respective objects according  SExtractor's segmentation map of the detection image.

 performs well on the two large spiral galaxies (IDs 19 and 22), recovering a close approximation of their shape, despite the symmetry and monotonicity constraints. 
While there will be cases where those constraints are too restrictive, e.g. for tightly wound spirals, these two examples demonstrate that galaxies with a modest amount of internal structure can be modeled well, even with a single component.
The combination of constraints is useful for dense groups of objects, e.g. around IDs 20 or 6, where e.g. a symmetry constraint alone would not suffice.
ID 6 is well modeled despite being rather faint because the simultaneous use of all available bands helps determine a best-fitting morphology, which in turn is used as a matched filter to determine the SED.
However, in IDs 20 or 25 we can see that this process can fail when undetected sources are present in the vicinity of model components.
The optimization algorithm of  will reduce the residuals by connecting to undetected emission via a flat “bridge”, which is the only option available under the monotonicity constraint.
If an undetected source has a different SED, a single-component source has no other option than to adopt the flux-weighted mean SED, which is clearly visible in ID 25.
Anℓ_0sparsity penalty (cf. <ref>) could suppress those bridges, however its use is problematic because it would bias low the recovery of total fluxes.



  §.§.§ Comparison to other deblenders
 SExtractor makes no such structural assumptions, but the limitations of the “cookie cutter” approach clearly show with the large spirals IDs 19 and 22, as well as their neighbors, e.g. ID 20.
As a result, most objects are too compact, which will lead to low biases on total fluxes and sizes.
 SExtractor is also prone to merging different sources together, even when they are discontinuous in the segmentation map (IDs 19, 22, 25).

The HSC-SDSS deblender performs surprisingly badly on ID 19 and, by extension, on its neighbors.
Because the pixel fluxes are exactly preserved, a failure in one source has strong effects on others, e.g. IDs 20 and 25 pick up the flux from the outskirts of 19.
However, the addition of monotonicity, which was not present in the SDSS  Photo implementation by Lupton (in prep.), helped in the three-in-a-row case of ID 6.

As a summary, we find  to perform well on a complex blended scene of deep HSC imaging, despite its simplifying assumptions of single-component sources and the rigidity of the symmetry and monotonicity constraints.
It simultaneously models all sources and is automatically consistent across all bands.
A conceptual limitation remains, namely the sensitivity to undetected sources or multiple stellar populations within detected sources, which could be mitigated by adding or increasing penalties such as a strongℓ_0sparsity that restricts the support of each component to a minimum number of pixels.
On the other hand, an iterative sequence that adds components at the location of significant and coherently colored residuals appears viable and could substantially augment traditional detection methods.



 §.§ Simulations


As the problem of deblending astronomical scenes is interrelated with detecting all relevant sources,
we seek to decouple these two tasks with the help of simulations where we know the number and location of each source.
The simulations also allow us to determine the accuracy of the recovered source properties, namely total fluxes, SEDs, and morphologies.
We will first test how well the assumptions of monotonicity and symmetry allow us to recover complex morphologies of isolated sources (<ref>) and then how that recovery is affected when multiple overlapping objects are present in each scene (<ref>).
Both sets of simulations are created with the  galsim software package <cit.>.

As metrics to quantify the fidelity of the source modeling and deblending, we determine the fluxes, SEDs, and per-pixel morphologies of deblender results.
For the total flux, which for  is given by∫d^2x  _k, we simply compute the fractional error compared to the input model.
For SEDs and morphologies we compute the correlation coefficient between the true and the observed vectorsv,

    ξ_v = v_true· v_obs/√(v_true· v_true)·√(v_obs· v_obs).

The coefficientsξ__kandξ__kmeasure how well a deblender recovers colors and shapes (including intensity) of componentk, respectively.
We consider these fidelity measures as a function of true source flux and source “blendedness”<cit.>
    β = 1-_k·_k/_all·_k,

where_all = ∑_k _k, expecting that good results will be harder to achieve for faint sources in very crowded areas.



  §.§.§ Individual galaxies

    < g r a p h i c s >
An example galaxy from the COSMOS  HST dataset in  galsim. 
    Top left: The “truth” image, defined by an original  HST image from which we removed the noise in the galaxy outskirts. 
    Top right: The difference between the original  HST image and the “truth”.
    Bottom left: The model generated by  with symmetry and monotonicity constraints.
    Bottom right: The residuals (truth - model) reveal spiral structure that is missing in the  model.     
    The panels on the left use a sinh^-1 scaling. 
    

We start by investigating how well  can reproduce realistic galaxy morphologies.
We use a set of galaxies from the COSMOS  HST sample described in <cit.> that are available with  galsim[<https://github.com/GalSim-developers/GalSim/wiki/RealGalaxy%20Data>].
From the full sample of 87,000 galaxies withi<25.2, we further select those with total fluxes between 300 and 2,000 counts that pass most stringent “marginal” cuts.
The bright cutoff was chose to remove galaxies with visible blending: many of the brightest galaxies have faint background galaxies or foreground stars in their wings.
As we seek to test the recovery of individual sources here, those interlopers are contaminants that reduce the correlation measure in <ref>.
They could easily be fit with  if their existence and positions were known, but for the simulations here we want to assume that we know the position of every individual source.
The lower flux cutoff was imposed to eliminate faint sources, whose outskirts are dominated by correlated noise from the original  HST observations.
To further reduce the effect of correlated noise, we estimate the noise level over the postage stamp and subtract the corresponding value from the entire image, setting all negative pixel values to zero.
While this procedure is likely to cut off some flux from the wings of the galaxy, the bulk of the galaxy and especially its internal structure is preserved.
The remaining sample comprises 799 galaxies.
An example galaxy from this procedure is shown in <ref>.


    < g r a p h i c s >
Morphology correlation coefficient (<ref>) of  models for 799 isolated galaxies from the COSMOS  HST sample in  galsim.

We run  on single-band images of all 799 galaxies. The pixel resolution is set to a typical ground-based value of 0.2 arcsec; no noise is added to the image.
We show the resulting model and the residuals for the example galaxy in the bottom panels of <ref>.
For reference, its total flux error is 0.08%, its morphology correlation coefficientξ_=0.999.
The RMS of the total flux error for the entire sample is 0.39%, which increases for fainter sources, suggesting that some correlated noise still made it into the input images and is suppressed by the monotonicity constraint.
The morphology correlationξ_is shown as a cumulative distribution function in <ref>.
These results demonstrate that for realistic galaxies,  recovers morphologies with a very strong correlation to the input morphologies despite the constraints of symmetry and monotonicity; consequently, the estimation of total fluxes is also excellent.



  §.§.§ Blended scenes


We now create a set of 10,000 blended scenes in the LSST filters  ugrizy, with widths and heights varying between 30 and 70 pixels, comprising on average a mix of 20% stars and 80% galaxies, with a common PSF across all bands.
The stellar sources have SEDs taken from the <cit.> templates for O5V, B5V, A5V, F5V, K5V, and M5V stars, and each stellar morphology is defined as a 2D Gaussian withσ=10^-9pix that is convolved with the PSF in each band.
Galaxy sources have SEDs based on four empirical <cit.> templates available in  galsim, randomly perturbed to generate a total of 25 different galaxy SED templates.
Their redshifts are randomly chosen between 0 and 2.1.
For the galaxy morphologies we use the same 799 galaxies selected in <ref>.
All sources are inserted into the scenes with a power-law distribution in radius from the center, resulting in strongly crowded blends in the middle of the images.
Uncorrelated Poisson noise is added to the images, using a sky background brightness at the mean noise level of the  HST images, so that we can largely hide the truncation of the input galaxy morphologies (cf. <ref>).


We run  and the current HSC-SDSS deblender and separate the results for stars and galaxies assuming that larger average brightness and more compact shape should make stars easier to deblend than galaxies.
Indeed, we find in <ref> that total flux errors are smaller for stars than for galaxies and that  substantially outperforms the HSC-SDSS deblender.
These two trends continue with the correlation coefficients of SEDs in <ref> and of morphologies in <ref>.
We understand the remarkable difference between these two deblending strategies, which have very similar constraints on the source morphologies, as a result of fitting a joint model across all bands () or independent models for each band (HSC-SDSS).
Our strategy is superior in these test cases because we combine images of all bands, resulting in a morphology that is driven by the band with the highest significance, which is in turn used as an optimal filter to determine the source SED. 

We emphasize that these results are achieved with single-component, symmetric, and monotonic source models, even though the true morphologies in the simulations exhibit more complex structure.
There will be cases in which those constraints lead to substantial errors, but we believe that they are now demonstrably justified as powerful default configurations for deblending in ground-based wide-field surveys.

We finally note that the runtime of  in our tests is about 10 ms per band and source, driven almost entirely by the cost to evaluate the likelihood gradients in <ref>.
The computational cost increases substantially if the PSF convolution cannot be expressed efficiently.
While slower than e.g. the HSC-SDSS deblender,  only has to be run once over all images and directly yields measurements of flux, SED, and morphology of the sources.
We fully expect source extraction with  to be a computationally expensive but still feasible task for deep optical surveys.


    < g r a p h i c s >
Comparison of the relative error of total fluxes in the g between  (top) and the HSC-SDSS deblender (bottom), for galaxies (left) and stars (right).
    For reference, the red horizontal lines show 5% relative errors.
    We only show the results of the r-band as the other bands look qualitatively similar.
    
    < g r a p h i c s >

    < g r a p h i c s >

Correlation, as defined by <ref>, between true SEDs and the SEDs from the HSC-SDSS deblender or   as a function of true source flux (top four panels) and blendedness (<ref>, bottom four panels.)
    < g r a p h i c s >

    < g r a p h i c s >

Same as <ref> for the correlation between the true source morphology and the morphologies from the HSC-SDSS deblender or . For HSC-SDSS, the morphology is taken from the g-band; other bands look similar.

 §.§ Extended Emission Line Regions of AGN

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

    < g r a p h i c s >

1st panel: False-color image of grizy coadd images of SDSSJ023106-034513 from the HSC Wide data release, shown with an arcsinh stretch. 2nd panel: Convolved  model of the scene with a two-component monotonic source for the host galaxy, a free-form source for the EELR, and monotonic and symmetric sources for other detections in the scene. 3rd panel: Partially deconvolved EELR component. 4th panel: Result of the band subtraction i-z from <cit.>, which subtracts the galactic continuum to isolates the EELR emission. 


To further demonstrate the capabilities of  as a flexible component separation framework, we analyze HSC Wide coadd images ingrizybands of galaxies with extended emission-line regions (EELRs) created by active galactic nuclei (AGN). 
Resulting from AGN radiation on inhomogeneous inter-stellar medium, the morphology of EELR is often complex, asymmetric, and non-monotonic. 
The strong emission lines, e.g., [OIII]λ5007, emitted can be detected and resolved by broadband imaging as components that may spatially overlap with but have colors distinct from their host galaxies <cit.>. 

The particular example shown in <ref>, SDSS J023106-034513, is taken from the EELR sample of <cit.>, who selected EELR candidates from spectroscopically identified type 2 obscured AGNs. 
With a high luminosity ofL_bol = 10^45.7erg s^-1, this AGN creates extended [OIII]λ5007 emission with a diameter of 60 kpc that is captured in the HSCi-band image, see <ref>. 

We seek to separate the EELR from its host galaxy with .
We first detect all sources in the image, and then assign them to either host, EELR, or one of 3 likely unrelated interlopers.
We constrain the morphology of the interlopers with our default combination of non-negativity, monotonicity and symmetry.
The host galaxy does not appear symmetric, so we remove that constraint[In experiments, we found that the host is better modeled by two co-centered components, similar to a bulge-disk model.],
while the EELR has no constraints other than non-negativity.

To compensate for the weakening of morphological constraints, we can utilize our knowledge of the EELR broadband SED based on the emission lines in the observed SDSS spectrum.
Similarly, by assuming that the host galaxy is responsible for the continuum spectrum without emission lines, we can also predict its broadband SED.
The proximal operator forA_kis thus simply a projection onto the precomputed SEDs of EELR or host and mute for any of the interlopers:

    _AGN(A_k) = A_EELR   if k = k_EELR
    A_Host   if k = k_Host
    A_k    else.

We use inverse per-pixel variances as weights and match the observed PSFs between bands by deconvolving to an effective resolution of 0.35 arcsec, using the formalism of <ref>.
We show the converged, convolved  model in the second panel of <ref> and the partially deconvolved EELR component in the third panel.
For comparison, in the forth panel, we also show the EELR image obtained via image continuum subtraction from <cit.>. 
In this approach, to isolate the EELR signal, the galaxy continuum in thei-band is removed by subtracting the (scaled)z-band image, which contains only the continuum but not line emissions. 
We can see that the  model of the EELR strongly resembles the continuum-subtraction result, while revealing smaller spatial features of the EELR.
A more detailed application of  to EELR extraction is planned for future work.



§ SUMMARY AND OUTLOOK


The employment of multi-channel information is long known to aid source separation schemes.
It was therefore expected that multi-band imaging should be beneficial in cases when multiple sources overlap in a celestial scene.
To exploit the availability of multi-band imaging from modern sky surveys, we have developed the source separation framework .
It uses a matrix factorization scheme, in which one matrix factor describes the amplitude of all components in each band and the other factor describes their spatial shape.
We use the proximal optimization method bSDMM presented by <cit.> to allow for an arbitrary number of constraints on each matrix factor.
The astrophysical interpretation of this model is that the scene is composed of a number of components that each describe one star or distinct stellar population.
Galaxies can be modeled as several such components to characterize e.g. bulge and disc-like features.

To aid the source separation further, especially in crowded regions with little variation in source colors (such as galaxy cluster cores), we constrain the component morphologies to be monotonic and/or symmetric with respect to the sources peak position.
If available, constraints on the SED of the sources, e.g. from stellar SED libraries or high-resolution spectra, are beneficial as well.

We have derived the analytical treatment of correlated pixel noise and of band-dependent PSF convolution, rendering  capable of working on coadd images with variable seeing conditions.
In essence, it performs a PSF-matched photometry measurement with an optimally chosen weight function given by the signal-to-noise-weighted mean morphology in all available bands.
Applied to an example scene from the HSC UltraDeep survey, we find   to robustly recover all detected sources, even if they exhibit complex morphologies.
For detailed analyses of extended galaxies one can conceive of additional constraints, such as the Fourier mode extension of  Galfit<cit.>.

We show with a set of dedicated simulations comprised of stars and realistic galaxy morphologies that  yields an accurate recovery of source fluxes, SEDs, and morphologies, even for prominent levels of blendedness.
In our tests it clearly outperforms the HSC-SDSS deblender in each of these measurements and is especially powerful for stars and marginally resolved as well as faint sources.

Because of its non-parametric nature,  exhibits a set of inherent advantages and drawbacks, which we summarize in <ref>. 
In addition, due to its flexibility it is sensitive to undetected sources or additional components of detected sources.
On the other hand, the inspection of the fit residuals appears to clearly reveal their presence.
We therefore envision a scheme in which components are added at the location of significant and coherently colored residuals.
A decision whether any newly added component is part of a previously detected source or independent of all other sources will require a post-processing step that encodes an interpretation based on the nature of the astronomical scene.

When run in this fashion  breaks a common model of data processing in optical astronomy that consists of three distinct stages: detection, deblending, and characterization.
If components are added from an analysis of the residuals and the models are used to infer source properties, it essentially covers aspects of all three stages.
While doing so, it makes assumptions about the SEDs and morphologies of the sources as well as the noise and the PSFs in the observations, which need to be valid to produce reliable results.
It will therefore be necessary to assess the accuracy of  models for specific science cases.
To facilitate inference directly on  sources, we intend to provide an error-estimation scheme, which accurately captures component degeneracies and effects of non-differentiable constraints.

We have intentionally developed  such that it can utilize any constraint or prior that can be expressed as a proximal operator without having to change the algorithm itself, thus rendering it suitable to a wide range of source separation problems in image data.
We demonstrated this flexibility with an AGN jet–host galaxy separation in this work and with a hyperspectral unmixing task in <cit.>.
It is inherent to the matrix factorization approach that  should work well in applications where high-quality information on the spatial distribution of sources is to be combined with high-quality information on their spectra.
An example is the combination of a single high-resolution image from a space-born telescope, such as  HST or the future  WFIRST<cit.>, with multi-band imaging of lower spatial resolution from a ground-based telescope.
We will pursue this research direction in forthcoming works.

 is implemented in  Python with  C++ extensions, accessible through  pybind11, and is available at <https://github.com/fred3m/scarlet>.
Its interface is highly modular and should be easily extendible to different applications.
Contributions and requests from the astronomical community are highly encouraged.

Main advantages and drawbacks of .
 Feature    Advantage   Disadvantage

Matrix Factorization     Low-rank representation of celestial scenes     Not suited for correlated changes of SED and morphology
0pt3exNon-parametric model     Arbitrary morphologies and SEDs     Requires critical sampling 
       No analytic PSF convolution
0pt3exProximal gradient optimization     Wide range of priors and penalties     No analytic error propagation



§ ACKNOWLEDGEMENTS

We gratefully acknowledge support from NASA grant 14-WPS14-0003.

This material is based upon work supported in part by the National Science Foundation through
Cooperative Agreement 1258333 managed by the Association of Universities for Research in Astronomy (AURA), and the Department of Energy under Contract No. DE-AC02-76SF00515 with the SLAC National
Accelerator Laboratory. Additional LSST funding comes from private donations, grants to universities,
and in-kind support from LSSTC Institutional Members.



§ REFERENCES
Miyazaki2018Komiyama2018Furusawa2018]