---
abstract: |
    Calcium imaging has revolutionized systems neuroscience, providing the ability to image large neural populations with single-cell resolution. The resulting datasets are quite large (with scales of TB/hour in some cases), which has presented a barrier to routine open sharing of this data, slowing progress in reproducible research. State of the art methods for analyzing this data are based on non-negative matrix factorization (NMF); these approaches solve a non-convex optimization problem, and are highly effective when good initializations are available, but can break down e.g. in low-SNR settings where common initialization approaches fail.

    Here we introduce an improved approach to compressing and denoising functional imaging data. The method is based on a spatially-localized penalized matrix decomposition (PMD) of the data to separate (low-dimensional) signal from (temporally-uncorrelated) noise. This approach can be applied in parallel on local spatial patches and is therefore highly scalable, does not impose non-negativity constraints or require stringent identifiability assumptions (leading to significantly more robust results compared to NMF), and estimates all parameters directly from the data, so no hand-tuning is required. We have applied the method to a wide range of functional imaging data (including one-photon, two-photon, three-photon, widefield, somatic, axonal, dendritic, calcium, and voltage imaging datasets): in all cases, we observe $\sim$2-4x increases in SNR and compression rates of 20-300x with minimal visible loss of signal, with no adjustment of hyperparameters; this in turn facilitates the process of demixing the observed activity into contributions from individual neurons. We focus on two challenging applications: dendritic calcium imaging data and voltage imaging data in the context of optogenetic stimulation. In both cases, we show that our new approach leads to faster and much more robust extraction of activity from the video data.
author:
- |
    E. Kelly Buchanan[^1],$^,$  Ian Kinsella,$^,$  Ding Zhou,$^,$  Rong Zhu[^2],   Pengcheng Zhou,\
    Felipe Gerhard[^3], John Ferrante,\
    Ying Ma[^4], Sharon H. Kim, Mohammed A Shaik,\
    Yajie Liang[^5], Rongwen Lu,\
    Jacob Reimer[^6], Paul G Fahey, Taliah N Muhammad,\
    Graham Dempsey, Elizabeth Hillman, Na Ji, Andreas S Tolias, Liam Paninski
bibliography:
- 'axon\_pipeline.bib'
title: 'Penalized matrix decomposition for denoising, compression, and improved demixing of functional imaging data'
---

Introduction {#introduction .unnumbered}
============

Functional imaging is a critical tool in neuroscience. For example, calcium imaging methods are used routinely in hundreds of labs, generating large-scale video datasets whose characteristics (cell shapes, signal-to-noise levels, background activity, signal timescales, etc.) can vary widely depending on the imaging modality and the details of the brain region and cell types being imaged. To handle this data, scientists must solve two basic tasks: we need to extract signals from the raw video data with minimal noise, and we need to store (and share) the data. A number of papers have focused on the first task [@mukamel2009automated; @maruyama2014detecting; @pnevmatikakis2016simultaneous; @pachitariu2016suite2p; @friedrich2017multi; @inan2017robust; @Reynolds2017; @petersen2017scalpel; @zhou2018efficient; @Mishne2018]; however, somewhat surprisingly, very little work has focused on the second task.

For both of these tasks, it is critical to denoise and compress the data as much as possible. Boosting the signal-to-noise ratio (SNR) is obviously important for detecting weak signals, performing single-trial analyses (where noise cannot be averaged over multiple trials), and for real-time experiments (where we may need to make decisions based on limited data - i.e., averaging over time is not an option). The benefits of compression are perhaps less obvious but are just as numerous: compression would facilitate much more widespread, routine open data sharing, enhancing reproducible neuroscience research. Compression will also be critical for in vivo imaging experiments in untethered animals, where data needs to be transmitted wirelessly, making data bandwidth a critical constraint. Finally, many signal extraction methods based on matrix factorization can be sped up significantly if run on suitably compressed data.

Previous methods for denoising and compressing functional data have several drawbacks. Generic video compression approaches do not take advantage of the special structure of functional imaging data and produce visible artifacts at high compression rates; more importantly, these approaches do not denoise the data, since they focus on compressing the full data, including noise, whereas our goal here is to discard the noise. Conversely, generic image denoising approaches do not offer any compression (and also fail to take advantage of strong structured correlations in the video data). Constrained nonnegative matrix factorization (CNMF) [@pnevmatikakis2016simultaneous] approaches provide state of the art denoising and demixing of calcium imaging data, but these methods can leave significant visible signal behind in the residual (discarding potentially valuable signal) and are highly dependent on the initialization of the matrix factorization; thus it would be dangerous to keep only the matrix factorization output and discard the raw data. Principal components analysis (PCA) is often employed as a compression and denoising method [@mukamel2009automated; @pachitariu2016suite2p], but PCA is based on a rather unstructured signal model and therefore provides a suboptimal encoder of functional data (we will discuss this point in further depth below). In addition, the computation time of PCA scales quadratically with the number of pixels (assuming a long video dataset) and therefore naive applications of PCA are rather slow [@friedrich2017multi]. Finally, importantly, it is difficult to automatically choose the number of principal components that should be retained in a given video (and the “correct” number of components can vary widely across different datasets).

Here we introduce a new simple approach to denoising and compressing functional video data. We apply a variant of penalized matrix decomposition [@Witten2009pmd] that operates locally in space, and encourages smoothness in both the spatial and temporal dimensions. This method offers multiple advantages over previous approaches. It is based on a signal model that is well-matched to the structure of the data: cells are local in space, there aren’t too many of them compared to the number of pixels (leading to a low-rank signal model), and cellular activity is smoother than the dominant noise sources, which are spatially and temporally uncorrelated. The approach is scalable (scaling linearly in the number of frames and pixels), and has modest memory requirements (because all processing is only performed in local spatial patches). All parameters (including the local matrix rank and the degree of smoothness of the output) are chosen automatically. Empirically we find that the method is highly effective, leaving behind minimal visible structure in the residual, while achieving 20-300x compression rates and 2-4x improvements in SNR. We demonstrate the method’s effectiveness on a wide variety of functional imaging datasets (both calcium and voltage imaging; one-, two- and three-photon imaging; and data including somas and dendrites) and show that the method is also effective on wide-field imaging data, where single-cell resolution is not available. Finally, we develop a new constrained NMF approach based on the denoised and compressed representation of the data, and apply this new demixing method to two challenging applications: dendritic calcium imaging data and voltage imaging data in the context of optogenetic stimulation. In both cases, we show that our new approach leads to faster and much more robust extraction of activity from the video data.

Methods {#methods .unnumbered}
=======

We begin by defining notation. Our starting point is an imaging dataset that has been motion-corrected (i.e., we assume that there is no motion of visible cellular components from frame to frame of the movie) and then “unfolded" into a $d \times T$ matrix $\mathbf{Y}$, where $T$ is the number of frames in the movie and $d$ is the number of pixels per frame (or voxels per frame if we are performing imaging in three dimensions). Now the typical approach is to model the data $\mathbf{Y}$ as $\mathbf{Y} = \mathbf{AC} + \mathbf{B} + \mathbf{E}$, where the columns of $\mathbf{A} \in \mathbb{R}^{d \times K}$ model the locations of each source (with $K$ sources total), the rows of $\mathbf{C} \in \mathbb{R}^{K \times T}$ model the time-varying fluorescence of each source, $\mathbf{B} \in \mathbb{R}^{d \times T}$ is a “background" term to handle signals that can not easily be split into single-neuronal components, and $\mathbf{E} \in \mathbb{R}^{d \times T}$ denotes temporally and spatially uncorrelated noise.

It is useful to break the processing pipeline into three sub-problems:

1.  **Denoising**: separation of neural signal $\mathbf{Y}^{*} = \mathbf{A}\mathbf{C} + \mathbf{B}$ from noise $\mathbf{E}$;

2.  **Compression** of signal $\mathbf{Y}^{*}$;

3.  **Demixing**: factorization of $\mathbf{Y}^{*}$ into its constituent components $\mathbf{A},\mathbf{C}$, and $\mathbf{B}$.

Most prior work has attempted to solve these sub-problems simultaneously, e.g., to recover $\mathbf{A}$ and $\mathbf{C}$ directly from the raw data $\mathbf{Y}$. As emphasized above, this direct approach involves a challenging non-convex optimization problem; the solution to this problem typically misses some structure in $\mathbf{Y}$, is highly sensitive to initialization and hyperparameter settings, and can be particularly unstable in low-SNR regimes. We have found empirically that a sequential approach is more robust and effective. First we compute the compressed and denoised estimate $\hat{\mathbf{Y}} = \mathbf{UV}$; here $\mathbf{U}$ and $\mathbf{V}$ are chosen so that $\hat{\mathbf{Y}}$ captures all of the signal in $\mathbf{Y}$ while retaining minimal noise (i.e., $\hat{\mathbf{Y}} \approx \mathbf{Y}^{*}$) and also $\mathbf{U}$ and $\mathbf{V}$ are highly-structured, compressible matrices, but we do not enforce any constraints between $(\mathbf{U}, \mathbf{V})$ and $(\mathbf{A}, \mathbf{C}, \mathbf{B})$. The computation of $\mathbf{U}$ and $\mathbf{V}$ essentially solves sub-problems 1 and 2 simultaneously. Second, we exploit $\mathbf{U}$, $\mathbf{V}$, and the resulting denoised $\hat{\mathbf{Y}}$ to facilitate the solution of problem 3. We discuss each of these steps in turn below.

Denoising & Compression {#denoising-compression .unnumbered}
-----------------------

To achieve good compression and denoising we need to take advantage of three key properties of functional imaging data:

1.  Signal sources are (mostly) spatially local;

2.  Signal is structured both temporally and spatially, whereas noise is temporally and spatially uncorrelated;

3.  Signal is (mostly) low-rank.

Given these structural assumptions, it is natural to construct $\mathbf{U}$ and $\mathbf{V}$ via a local penalized matrix decomposition approach[^7]: we break the original data matrix $\mathbf{Y}$ into a collection of overlapping spatial patches, then decompose each of these matrix patches (in parallel) using a factorization method that enforces smoothness in the estimated spatial and temporal factors, then combine the resulting collection of spatial and temporal factors over all the patches into a final estimate of $\mathbf{U}$ and $\mathbf{V}$. (See [CaImAn](https://github.com/flatironinstitute/CaImAn) for a similar patch-wise approach to the demixing problem.)

We have experimented with several approaches to penalized matrix decomposition (PMD), and found that an iterative rank-one deflation approach similar to the method described in [@Witten2009pmd] works well. We begin by standardizing the data within a patch: for each pixel, we subtract the mean and normalize by an estimate of the noise variance within each pixel; the noise variance is estimated using the frequency-domain method described in [@pnevmatikakis2016simultaneous], which exploits the fact that the signal and noise power are spectrally separated in movies with sufficiently high frame rates. After this normalization we can model the noise $\mathbf{E}$ as roughly spatially and temporally homogeneous. Denote this standardized data matrix within a patch as $\mathbf{Y_0}$, and Frobenius norm as $||.||_F$. Then at the $k^{th}$ iteration PMD extracts the best rank-one approximation $\mathbf{u}_k\mathbf{v}_k^T$ to the current residual $\mathbf{R}_k = \mathbf{Y_0} - \sum_{n=1}^{k-1} \mathbf{u}_n\mathbf{v}_n^T$, as determined by the objective $$(\mathbf{u}_k, \mathbf{v}_k) = \underset{\mathbf{u}, \mathbf{v}}{\arg\min} ~ || \mathbf{R}_k - \mathbf{u}  \mathbf{v}^T ||_F \hspace{1em} \text{subject to} \hspace{1em}
P_{spatial}(\mathbf{u}) \leq c_{1}^k,\  P_{temporal}(\mathbf{v}) \leq c_{2}^{k}, \label{eqn:CPMD}$$ followed by a temporal debiasing update $\mathbf{v}_k = \mathbf{R}_k^T\mathbf{u}_k$. The objective (\[eqn:CPMD\]) can be ascended via alternating minimization on $\mathbf{u_k}$ and $\mathbf{v_k}$.

Note that if we drop the $P_{spatial}(\mathbf{u})$ and $P_{temporal}(\mathbf{v})$ constraints above then we can solve for $\mathbf{u}_k$ and $\mathbf{v}_k$ directly by computing the rank-1 singular value decomposition (SVD) of $\mathbf{R}_k$; in other words, by performing PCA within the patch. Since we have normalized the noise scale within each pixel, PCA should identify the signal subspace within the patch, given enough data (because the normalized projected data variance in any direction will be equal to one plus the signal variance in this direction; since PCA searches for signal directions that maximize variance, PCA will choose exactly the signal subspace in the limit of infinite data). Indeed, as discussed in the results section, simple patch-wise PCA (with an appropriate adaptive method for choosing the rank) often performs well, but incorporating spatial and temporal penalties in the optimization can push $\mathbf{u}_k$ and $\mathbf{v}_k$ closer to the signal subspace, resulting in improved compression and SNR.

How should we define the penalties $P_{spatial}(\mathbf{u})$ and $P_{temporal}(\mathbf{v})$, along with the corresponding constraints $c_{1}^k$ and $c_{2}^{k}$? The simplest option would be to use quadratic smoothing penalties; this would lead to a simple closed-form linear smoothing update for each $\mathbf{u_k}$ and $\mathbf{v_k}$. However, the signals of interest here have inhomogeneous smoothness levels — an apical dendrite might be spatially smooth in the apical direction but highly non-smooth in the orthogonal direction, and similarly a calcium signal might be very smooth except at the times at which a spike occurs. Therefore simple linear smoothing is typically highly suboptimal, often resulting in both undersmoothing and oversmoothing in different signal regions. We have found total variation (TV) [@Rudin:1992:NTV:142273.142312] and trend filtering (TF) [@Kim2009tf] penalties to be much more empirically effective. We let $$\begin{aligned}
P_{temporal}(\mathbf{v}) = \| \mathbf{D}^{(2)} \mathbf{v}\|_1 = \sum_{t=2}^{T-1} |\mathbf{v}_{t-1} - 2 \mathbf{v}_{t} + \mathbf{v}_{t+1}| \end{aligned}$$ and $$\begin{aligned}
P_{spatial}(\mathbf{u}) = 
\|\mathbf{\nabla}_{\mathcal{G}}\mathbf{u}\|_1 = \sum_{(i,j) \in \mathcal{E}} | \mathbf{u}_i - \mathbf{u}_j |.\end{aligned}$$ Here $\mathbf{D}^{(2)}$ denotes the one-dimensional discrete second order difference operator and $\mathbf{\nabla}_{\mathcal{G}}$ the incidence matrix of the nearest-neighbor pixel-adjacency graph (pixels $(i,j)$ are in the edge set $\mathcal{E}$ if the pixels are nearest neighbors).

Similarly to [@pnevmatikakis2016simultaneous], we define the smoothing constraints $c_1^k$ and $c_2^k$ implicitly within the alternating updates by the simple reformulation $$\mathbf{u}_k = \underset{\mathbf{u}}{\arg\min} \| \mathbf{R}_k \mathbf{v}_k - \mathbf{u}\|_2^2\  s.t.\  \| \mathbf{\nabla}_{\mathcal{G}} \mathbf{u}\|_1 \leq c_1^k
\iff 
\mathbf{u}_k = \underset{\mathbf{u}}{\arg\min} \| \mathbf{\nabla}_{\mathcal{G}} \mathbf{u}\|_1\  s.t.\ \| \mathbf{R}_k \mathbf{v}_k - \mathbf{u}\|_2^2  \leq \hat{\sigma}^2_{\tilde{\mathbf{u}}} d \label{eqn:spatial_update}$$ and $$\mathbf{v}_k = \underset{\mathbf{v}}{\arg \min} \| \mathbf{R}_k ^T \mathbf{u}_k - \mathbf{v}\|_2^2 \  s.t.\ \| \mathbf{D}^{(2)} \mathbf{v}\|_1 \leq c_2^k
\iff 
\mathbf{v}_k = \underset{\mathbf{v}}{\arg\min} \| \mathbf{D}^{(2)} \mathbf{v}\|_1 \ s.t.\ \| \mathbf{R}_k ^T \mathbf{u}_k - \mathbf{v}\|_2^2  \leq \hat{\sigma}^2_{\tilde{\mathbf{v}}} T \label{eqn:temporal_update}$$ where $\hat{\sigma}^2_{\tilde{\mathbf{u}}}$ (resp. $\hat{\sigma}^2_{\tilde{\mathbf{v}}}$) estimates the noise level of the unregularized update $\tilde{\mathbf{u}}_k = \mathbf{R}_k \mathbf{v}_k$ (resp. $\tilde{\mathbf{v}}_k = \mathbf{R}_k^T \mathbf{u}_k$), and we are using the fact that if the residual $\mathbf{R}_k \mathbf{v}_k - \mathbf{u}$ contains just noise then its squared norm should be close to $\hat{\sigma}^2_{\tilde{\mathbf{u}}} d$, by the law of large numbers (and similarly for equation \[eqn:temporal\_update\]). See Algorithm \[alg:ROD\] for a summary.

To solve the constrained problems on the right-hand side we use the line search approach described in [@Langer2017cps]. We solve the primal form of the TV optimization problem (\[eqn:spatial\_update\]) using the proxTV package [@barberoTV14], and of the TF optimization problem (\[eqn:temporal\_update\]) using the Primal-Dual Active Set method in [@Han2016pdas]. Both of these methods can exploit warm starts, leading to major speedups after a good initial estimate is found. Empirically the TF optimization scales linearly with the movie length $T$; since the scale of the TV problem is bounded (because we work in local spatial patches) we have not explored the scaling of the TV problem in depth.

Figure \[fig:TrendFiltering\] illustrates the effect of trend filtering on a couple $\mathbf{v}$ components. One important difference compared to previous denoising approaches [@haeffele2014structured; @pnevmatikakis2016simultaneous] is that the TF model is more flexible than the sparse autoregressive model that is typically used to denoise calcium imaging data: the TF model does not require the estimation of any sparsity penalties or autoregressive coefficients, and can handle a mixture of positive and negative fluctuations, while the sparse nonnegative autoregressive model can not (by construction). This is important in this context since each component in $\mathbf{V}$ can include multiple cellular components (potentially with different timescales), mixed with both negative and positive weights.

![Illustration of trend filtering. Each row shows a component $\mathbf{v}$ extracted from the voltage imaging dataset (see Results section for details). Red indicates simple projected signal $\tilde{\mathbf{v}} = \mathbf{R}^T \mathbf{u}$; blue indicates $\mathbf{v}$ after trend filtering. Errorbars on left indicate $2 \times$ estimated noise scale; right panels show zoomed region indicated by dashed lines in left panel.[]{data-label="fig:TrendFiltering"}](./methods/Effect_Of_TF.png){width="17cm"}

To complete the description of the algorithm on a single patch we need an initialization and a stopping criterion to adaptively choose the rank of $\mathbf{U}$ and $\mathbf{V}$. For the latter, the basic idea is that we want to stop adding components $k$ as soon as the residual looks like uncorrelated noise. To make this precise, we define a pair of spatial and temporal “roughness" test statistics $$\begin{aligned}
&T_{temporal}(\mathbf{v}) = \|\mathbf{D}^{(2)} \mathbf{v}\|_1 / \| \mathbf{v} \|_1
&T_{spatial}(\mathbf{u}) = \|\mathbf{\nabla}_{\mathcal{G}}\mathbf{u}\|_1 / \| \mathbf{u} \|_1\end{aligned}$$ and compute these statistics on each extracted $\mathbf{u}_k$ and $\mathbf{v}_k$. We accept or reject each component according to a one-sided hypothesis test under the null hypothesis that $\mathbf{R}_k$ consists of uncorrelated Gaussian noise of variance one. (We compute the critical region for this test numerically.) In the compression stage we are aiming to be rather conservative (we are willing to accept a bit of extra noise or a slightly higher-rank $\mathbf{U}$ and $\mathbf{V}$ in order to ensure that we are capturing the large majority of the signal), so we terminate the outer loop (i.e., stop adding more components $k$) after we reject a couple components $k$ in a row. See Algorithm \[alg-full-PMD\] for a summary.

To initialize, we have found that setting $\mathbf{u}_0 \propto \mathbf{1}$ works well. To speed up early iterations, it is natural to iterate the projections while skipping the denoising steps; this corresponds to intializing with an approximate rank-1 SVD as computed by power iterations. Initializing in this manner can reduce the total number of iterations needed for $\mathbf{u}_k, \mathbf{v}_k$ to converge. Matrix-vector multiplications are a rate limiting step here; thus, these initial iterations can be sped up using spatial and temporal decimation on $\mathbf{R}_k$. Empirically, decimation has the added benefit of boosting signal (by averaging out noise in neighboring timepoints and pixels) and can be useful for extracting weak components in low SNR regimes; see [@friedrich2017multi] for a related discussion.

The method described so far handles a single spatial patch of data. We can process patches in parallel; a multi-core implementation of this method (assigning different patches to different cores) achieves nearly linear speedups. We have found that for some datasets edge artifacts can appear near patch boundaries if the patches do not overlap spatially. These boundary artifacts can be eliminated by performing a $4 \times$ over-complete block-wise decomposition of $\mathbf{Y}$ using half-offset grids for the partitions (so that each pixel $x$ lies within the interior of at least one patch). Then we combine the overlapping patches together via linear interpolation (see [@pnevmatikakis2017normcorre] for a similar approach): set $$\hat{\mathbf{Y}}(x,t) = \frac{\sum_{p} \mathbf{a}_p(x) \hat{\mathbf{Y}}_p(x,t)} {\sum_p \mathbf{a}_p(x)},$$ where $p$ indexes the patches (so $\hat{\mathbf{Y}}_p$ denotes the denoiser output in the $p$-th patch) and $0 \leq \mathbf{a}_p(x) \leq 1$ is a “pyramid" function composed of piecewise linear functions that start at $0$ at the patch boundaries and increase linearly to $1$ at the center of the patch.

The above is equivalent to starting with a collection of overlapping sparse local factorizations $\mathbf{U}_p \mathbf{V}_p$, forming element-wise products between the individual spatial components $\mathbf{U}_{ip}$ and the pyramid functions $\mathbf{a}_p$, and then forming the union of the result to obtain a new factorization $\mathbf{UV}$. Typically this will result in some redundancy due to the overlapping spatial components; we remove this redundancy in a final backwards model selection step that tests whether each temporal component can be explained as a weighted sum of its neighbors. More precisely, we sort the components in ascending order according to the $L_2$ norms of $\mathbf{U}_{ip} \cdot a_p$. For each $i$ in this order we then regress $\mathbf{V}_i$ onto the collection of temporal components $\mathbf{V}_j$ whose corresponding spatial components $\mathbf{U}_j$ overlap with $\mathbf{U}_i$, i.e., approximate $ \hat{\mathbf{V}_i} = \sum_j \beta_j \mathbf{V}_j$. We then test the signal strength of the residual $\mathbf{V}_i - \hat{\mathbf{V}_i}$ (using the temporal test statistic defined previously); the component is rejected if the residual is indistinguishable from noise according to this test statistic. If component $i$ is rejected then we distribute its energy to the remaining spatial components according to the regression weights: $\mathbf{U}_{j} = \mathbf{U}_{j} + \beta_{j} \mathbf{U}_{i}$.

We conclude with a few final implementation notes. First, the results do not depend strongly on the precise patch size, as long as the patch size is comparable to the spatial correlation scale of the data: if the patches are chosen to be much smaller than this then the $\mathbf{V}$ components in neighboring patches are highly correlated, leading to excessive redundancy and suboptimal compression. (Conversely, if the patch size is too big then the sparsity of $\mathbf{U}$ is reduced, and we lose the benefits of patch-wise processing.)

Second, in some datasets (e.g., widefield imaging, or microendoscopic imaging data), large background signals are present across large portions of the field of view. These background signals can be highly correlated across multiple spatial patches, leading to a suboptimal compression of the data if we use the simple independent-patch approach detailed above. Thus in some cases it is preferable to run a couple iterations of PMD(TV, TF) on the full $\mathbf{Y}$ and then subtract the resulting components away before moving on to the independent block processing scheme. We have found that this effectively subtracts away dominant background signals; these can then be encoded as a small number of dense columns in the matrix $\mathbf{U}$, to be followed by a larger number of sparse columns (corresponding to the small patches), resulting in an overall improvement in the compression rate. See the for an example.

The patch-wise PMD(TV,TF) approach results in an algorithm that scales linearly in three critical parameters: $T$ (due to the sparse nature of the second-difference operator in the TF step), $d$ (due to the patch-wise approach), and the rank of $\mathbf{U}$ and $\mathbf{V}$. We obtain further speedups by exploiting warm starts and parallel processing over patches. Additional speedups can be obtained for very long datasets by computing $\mathbf{U}$ on a subset of the data and then updating $\mathbf{V}$ on the remainder of the movie; the latter step does not require any PMD iterations (since the spatial signal subspace has already been identified) and is therefore very fast, just requiring a single temporal update call per element of $\mathbf{V}$.

Demixing {#demixing .unnumbered}
--------

The methods described above provide a compressed and denoised representation of the original data $\mathbf{Y}$: the output matrices $\mathbf{U}$ and $\mathbf{V}$ are low-rank compared to $\mathbf{Y}$, and $\mathbf{U}$ is additionally highly sparse (since $\mathbf{U}$ is formed by appending spatial components $\mathbf{u}$ from multiple local spatial patches, and each $\mathbf{u}_k$ is zero outside of its corresponding patch). How can we exploit this representation to improve the demixing step?

It is useful to first take a step back to consider the strengths and weaknesses of current state of the art demixing methods, most of which are based on NMF. The NMF model is very natural in calcium imaging applications, since each neuron has a shape that is fixed over the timescale of a typical imaging experiment (and these shapes can be represented as non-negative images, i.e., an element of the $\mathbf{A}$ matrix), and a corresponding time-varying calcium concentration that can be represented as a non-negative vector (an element of $\mathbf{C}$): to form a movie we simply take a product of each of these terms and add them together with noise and background, i.e., form $\mathbf{Y}= \mathbf{AC} + \mathbf{B} + \mathbf{E}$.

However, current NMF-based approaches leave room for improvement in several key directions. First, since NMF is a non-convex problem, good initializations are critical to obtain good results via the standard alternating optimization approaches (similar points are made in [@petersen2017scalpel]). Good initialization approaches have been developed for somatic or nuclear calcium imaging, where simple Gaussian shape models are useful crude approximations to the elements of $\mathbf{A}$ [@pnevmatikakis2016simultaneous], but these approaches do not apply to dendritic or axonal imaging. Second (related), it can be hard to separate weak components from noise using current NMF-based approaches. Finally, voltage imaging data does not neatly fit in the NMF framework, since voltage traces typically display both positive and negative fluctuations around the baseline resting potential.

To improve the robustness of NMF approaches for demixing functional data, we make use of the growing literature on “guaranteed NMF” approaches — methods for computing a non-negative matrix factorization that are guaranteed to output the “correct” answer under suitable conditions and assumptions [@donoho2004does; @recht2012factoring; @arora2012computing; @li2016recovery]. In practice, these methods work well on clean data of sufficiently small dimensionality, but are not robust to noise and scale poorly to high-dimensional data. We can solve both of these issues by “superpixelizing" the denoised version of $\mathbf{Y}$; the resulting NMF initialization method improves significantly on state of the art methods for processing dendritic and axonal data. We also take advantage of the sparse, low-rank structure of $\mathbf{U}$ and $\mathbf{V}$ to speed up the NMF iterations.

### Initialization via pure superpixels {#initialization-via-pure-superpixels .unnumbered}

The first step of the initialization procedure is to identify groups of highly correlated spatially connected pixels – “superpixels." The idea is that a pixel within a neuron should be highly correlated with its neighbors, while a pixel containing mostly noise should have a much lower neighbor correlation. These neighbor correlations, in turn, can be estimated much more accurately from the denoised compared to the raw data. The superpixelization procedure results in a set of non-overlapping groups of pixels which are likely to be contained in good neural components. Then we want to extract “pure” superpixels, i.e., the subset of superpixels dominated by signal from just one neural component. We will use the temporal signals extracted from these pure superpixels to seed $\mathbf{C}$ in the NMF decomposition.

To identify superpixels, we begin with the denoised data $\hat{\mathbf{Y}} = \mathbf{UV}$. Since the compression process discussed in the previous section is rather conservative (aiming to preserve the full signal, at the expense of retaining a modest amount of noise), there is room to apply a more aggressive lossy denoiser in the initialization stage to further reduce any remaining noise in $\hat{\mathbf{Y}}$. We soft-threshold signals in each pixel that are not sufficiently large — less than the median plus $\delta \times$ the median absolute deviation (MAD) within each pixel, with $\delta \approx 1$ or $2$. (This thresholding serves to extract mostly spiking activity from functional imaging data.) We identify two neighboring pixels to be from the same superpixel if their resulting denoised, soft-thresholded temporal signals have a correlation larger than a threshold $\epsilon$, with $\epsilon \approx 0.9$. Superpixels that contain fewer than $\tau$ pixels are discarded to further reduce noise and the total number of superpixels. We then apply rank 1 NMF on the signals from each superpixel to extract their (thresholded) temporal activities.

To extract pure superpixels, we apply the Successive Projection Algorithm (SPA) [@gillis2014fast] to the temporal activities of superpixels. This algorithm removes “mixed” superpixels whose temporal activity can be modeled as a nonnegative linear combination of activity in other superpixels (up to some R-squared level larger than $ 1-\kappa$, where we use $\kappa \approx 0.2$) and outputs the remaining “pure" superpixels. See Algorithm \[alg1\] for pseudocode.

Note that running SPA on superpixels rather than raw pixels improves performance significantly here, since averaging signals within superpixels boosts SNR (making it easier to separate signal from noise and isolate pure from mixed pixels) and also greatly reduces the dimensionality of the non-negative regression problem SPA has to solve at each iteration. (To keep the problem size small we also run SPA just on small local spatial patches, as in the previous section.) Finally, while we have obtained good results with SPA, other approaches are available [@gillis2018fast] and could be worth further exploration in the future. See Figure \[fig:vi\_superpixels\] for a visual summary of the full procedure.

  -- --
     
     
     
     
     
     
     
     
  -- --

### Local NMF {#local-nmf .unnumbered}

Next we run NMF, using the temporal signals extracted from the “pure” superpixels to initialize $\mathbf{C}$. Given the initial $\mathbf{C}$, the typical next step is to regress onto the data to initialize $\mathbf{A}$. (Note that pure superpixels typically capture just a subset of pixels within the corresponding neuron, so it is not efficient to initialize $\mathbf{A}$ with the pure superpixels.) However, given the large number of pixels in a typical functional imaging video, direct regression of $\mathbf{C}$ onto $\mathbf{Y}$ is slow and overfits, providing poor estimates of $\mathbf{A}$.

This issue is well-understood [@pnevmatikakis2016simultaneous], and several potential solutions have been proposed. For somatic imaging it makes sense to restrict the support of $\mathbf{A}$ to remain close to their initial values (we could use a dilation of the superpixel support for this). But for data with large dendritic or axonal components this approach would cut off large fractions of these components. Sparse regression updates are an option here, but these do not enforce spatial structure in the resulting $\mathbf{A}$ directly; this often results in “speckle" noise in the estimated spatial components (c.f. Figure \[realcompare\] below).

We have found the following approach to be more effective. We initialize the support set $\Omega_k$ as the support of the $k$-th “pure” superpixel. Given $\mathbf{C}$, we compute the correlation image for each component $k$ as the correlation between the denoised data $\hat{\mathbf{Y}}$ and the $k$-th temporal component, $\mathbf{C}_k$. We truncate this correlation image below a certain threshold $\epsilon_1$ to zero, then update $\Omega_k$ as the connected component of the truncated correlation image which overlaps spatially with the previous $\Omega_k$. We use the modified fastHALS algorithm in [@friedrich2017multi] to update $\mathbf{A}$, $\mathbf{C}$, and $\mathbf{B}$ to locally optimize the objective $$\label{objnmf}
    \min_{\mathbf{A},\mathbf{C},{\textbf{\textit{b}}}} \|\hat{\mathbf{Y}}-\mathbf{AC} -\mathbf{B}\|_{F}^2,\ \mathrm{s.t.}\ \mathbf{A}_k^x = 0 ~ \forall x \not \in \Omega_k , \mathbf{A}\geqslant 0, \mathbf{C}\geqslant 0, \mathbf{B}={\textbf{\textit{b}}}\mathbf{1}^T, {\textbf{\textit{b}}}\geqslant 0.$$ Here we have modeled the background $\mathbf{B}$ as a simple temporally-constant vector; we discuss generalizations to time-varying backgrounds below. Also note that we are approximating $\hat{\mathbf{Y}}$ directly here, not the thresholded version we used to extract the superpixels above.

Finally, we incorporate a merge step: we truncate the correlation image below certain threshold $\epsilon_2$ to zero, and automatically merge neurons if their truncated correlation images are highly overlapped. The full algorithm is shown in Algorithm \[alg2\].

### Further implementation details {#further-implementation-details .unnumbered}

*Multi pass strategy:* As in [@zhou2018efficient], we find it effective to take a couple passes over the data; particularly in datasets with high neuron density, the first NMF pass might miss some dim neurons. We decrease the MAD threshold $\delta$ and re-run Algorithm \[alg1\] on the residual to find additional components, and then run a final merge and NMF update to complete the pipeline.

*Improvements from denoising and compression:* Compressed data leads to faster NMF updates, since we can replace $\hat{\mathbf{Y}}$ as $\mathbf{UV}$; in fastHALS, we can regress each ${\textbf{\textit{a}}}_k$ on $\mathbf{U}$ or ${\textbf{\textit{c}}}_{k}$ on $\mathbf{V}$ first instead of directly onto $\mathbf{Y}$. Similarly, when calculating the correlation image, we can compute the correlation between the low rank $\mathbf{V}$ and ${\textbf{\textit{c}}}_k$ first. As emphasized above, denoising also improves the estimation of the correlation images, which in turn improves the estimation of the support sets $\Omega_k$.

*Time-varying background:* It is straightforward to generalize the objective \[objnmf\] to include a time-varying background, using either a low-rank model (as in [@pnevmatikakis2016simultaneous]) or a ring-structured model (as in [@zhou2018efficient]). For the low-rank background model, we have found that performing an SVD on the data excluding the support of the superpixels provides an efficient initialization for the background temporal components.

*Incorporating temporal penalties*: Note that we are only imposing nonnegativity in $\mathbf{C}$ here; after denoising to obtain $\hat{\mathbf{Y}}$, we have found that this simple nonnegative constraint is sufficient for the datasets examined here. However, it is certainly possible to incorporate temporal penalties or constraints on $\mathbf{C}$ (e.g., a TF penalty or a non-negative auto-regressive penalty as in [@pnevmatikakis2016simultaneous]), either within each iteration or as a final denoising step.

*Post-processing*: We find that sorting the extracted components by their “brightness," computed as $\max {\textbf{\textit{a}}}_k\cdot\max{\textbf{\textit{c}}}_k$, serves to separate dim background components from bright single-neuronal components. We also found it useful to drop components whose temporal trace has skewness less than 0.5; traces with high skewness correspond to components with significant spiking activity, but low-skewness traces corresponded to noise.

Motion corrected data $\mathbf{Y}\in \mathbb{R}^{d\times T}$, MAD threshold $\delta$, minimum size of superpixels $\tau$, correlation threshold for superpixels $\epsilon$, $R^2$ threshold in SPA $\kappa$. $\sigma({\textbf{\textit{x}}})\leftarrow$ estimated noise for each pixel ${\textbf{\textit{x}}}$ of $\mathbf{Y}$; $\mu({\textbf{\textit{x}}})\leftarrow$ mean for each pixel of $\mathbf{Y}$; $\mathbf{Y} \leftarrow \left(\mathbf{Y}-\mu({\textbf{\textit{x}}})\right) / \sigma({\textbf{\textit{x}}})$; $(\hat{\mathbf{Y}},\mathbf{U},\mathbf{V}) \leftarrow$ PMD($\mathbf{Y}$); $n \leftarrow 0$; $\mathbf{A} \leftarrow [\ ]$, $\mathbf{C}\leftarrow [\ ]$, ${\textbf{\textit{b}}}\leftarrow\mathrm{median}$ for each pixel of $\hat{\mathbf{Y}}$; $\mathbf{R} \leftarrow \hat{\mathbf{Y}} -\mathbf{AC} - {\textbf{\textit{b}}}$; $\sigma_{med}({\textbf{\textit{x}}})\leftarrow$ median absolute deviation for each pixel of $\mathbf{R}$; $\mu_{med}({\textbf{\textit{x}}})\leftarrow$ median for each pixel of $\mathbf{R}$; $\tilde{\mathbf{Y}}\leftarrow \max\left(0, \mathbf{R} - \mu_{med}({\textbf{\textit{x}}}) - \delta\cdot\sigma_{med}({\textbf{\textit{x}}})\right)$; $\mathrm{corr}({\textbf{\textit{x}}},{\textbf{\textit{x}}}^*)\leftarrow \mathrm{corr}\left(\tilde{\mathbf{Y}}({\textbf{\textit{x}}},t),\tilde{\mathbf{Y}}({\textbf{\textit{x}}}^*,t)\right)$ for all neighbouring pixel pairs $({\textbf{\textit{x}}},{\textbf{\textit{x}}}^*)$; Extract superpixels: connect ${\textbf{\textit{x}}}$ and ${\textbf{\textit{x}}}^*$ together if $\mathrm{corr}({\textbf{\textit{x}}},{\textbf{\textit{x}}}^*)\geqslant \epsilon$ to construct connected components and discard those smaller than $\tau$, forming superpixels $\Omega_k,k=1,\cdots,K$; $({\textbf{\textit{a}}}_{k}, {\textbf{\textit{c}}}_{k})\leftarrow \mathrm{rank\ 1\ NMF}$ of $\tilde{\mathbf{Y}}$ on support $\Omega_k , k= 1,\cdots, K$; $[i_1,i_2,\cdots,i_S]\leftarrow\mathrm{SPA}([{\textbf{\textit{c}}}_{1},{\textbf{\textit{c}}}_{2},\cdots,{\textbf{\textit{c}}}_{K}], \kappa)$; $i_1,i_2,\cdots,i_S$ are indices of pure superpixels; $\mathbf{A}_0\leftarrow[\mathbf{A}, {\textbf{\textit{a}}}_{i_1},{\textbf{\textit{a}}}_{i_2},\cdots,{\textbf{\textit{a}}}_{i_S}]$; $\mathbf{C}_0\leftarrow[\mathbf{C}^T, {\textbf{\textit{c}}}_{i_1},{\textbf{\textit{c}}}_{i_2},\cdots,{\textbf{\textit{c}}}_{i_S}]^T$; ${\textbf{\textit{b}}}_0\leftarrow {\textbf{\textit{b}}}$; $(\mathbf{A}, \mathbf{C}, {\textbf{\textit{b}}})\leftarrow\mathrm{LocalNMF}(\mathbf{U}, \mathbf{V}, \mathbf{A}_0, \mathbf{C}_0, {\textbf{\textit{b}}}_0)$; $\delta \leftarrow \delta-1$; $n \leftarrow n+1$; $\eta(k)\leftarrow$ estimated noise for ${\textbf{\textit{c}}}_k$ using average of high frequency domain of PSD; (Optional) Denoise temporal components, e.g. by $\ell_1$ trend filter: ${\textbf{\textit{c}}}_k\leftarrow \min\limits_{\tilde{{\textbf{\textit{c}}}}_k} \|\tilde{{\textbf{\textit{c}}}}_k\|_1,\ \mathrm{s.t.}\ \|\tilde{{\textbf{\textit{c}}}}_k-{\textbf{\textit{c}}}_k\|_{F}\leqslant \eta(k)\sqrt{T}, k=1,\cdots,K$; $\mathbf{A},\mathbf{C},{\textbf{\textit{b}}}$

Compressed factors $\mathbf{U} \in \mathbb{R}^{d\times r}, \mathbf{V} \in \mathbb{R}^{T\times r}$ ($r = rank (\hat{\mathbf{Y}})$); initial constant background ${\textbf{\textit{b}}}_0$, spatial components $\mathbf{A}_0=[{\textbf{\textit{a}}}_{1,0},\cdots,{\textbf{\textit{a}}}_{K,0}]\in\mathbb{R}^{d\times K}$, and temporal components $\mathbf{C}_0=[{\textbf{\textit{c}}}_{1,0},\cdots,{\textbf{\textit{c}}}_{K,0}]^T \in\mathbb{R}^{K\times T}$; truncation threshold when updating support $\epsilon_1$, truncation threshold when merging $\epsilon_2$, overlap threshold when merging $\epsilon_3$. $\Omega_k \leftarrow \mathrm{supp}({\textbf{\textit{a}}}_{k,0})$ is spatial support for $k$-th component, $k=1,\cdots,K$; $\hat{\mathbf{A}} \leftarrow \mathbf{A}_0, \hat{\mathbf{C}}\leftarrow \mathbf{C}_0, \hat{{\textbf{\textit{b}}}}\leftarrow{\textbf{\textit{b}}}_0$; $\nu({\textbf{\textit{x}}})\leftarrow$ standard deviation for each pixel of $\hat{\mathbf{Y}} = \mathbf{UV}$; $\bar{\mathbf{V}}\leftarrow$ mean for each column of $\mathbf{V}$; $\mathbf{P} \leftarrow \left[\mathbf{U},-{\textbf{\textit{b}}}\right]\left(
        \begin{bmatrix}
        \mathbf{V}\\
        \mathbf{1}^T\\
        \end{bmatrix}\hat{\mathbf{C}}^T\right)$; $\mathbf{Q} \leftarrow \hat{\mathbf{C}}\hat{\mathbf{C}}^{T}$; Update spatial: $\hat{{\textbf{\textit{a}}}}_{k}(\Omega_k) \leftarrow \max\left(0, \hat{{\textbf{\textit{a}}}}_{k}(\Omega_k) + \frac{\mathbf{P}(\Omega_k,k)-\hat{\mathbf{A}}(\Omega_k)\mathbf{Q}(:,k)}{\mathbf{Q}(k,k)}\right)$; Update constant background: $\hat{{\textbf{\textit{b}}}} \leftarrow \max\left(0, \frac{1}{T}(\mathbf{UV}-\hat{\mathbf{A}}\hat{\mathbf{C}})\mathbf{1}\right)$; $\mathbf{P} \leftarrow \left[\mathbf{V}^T,\mathbf{1}\right]\left(\left[\mathbf{U},-{\textbf{\textit{b}}}\right]^T\hat{\mathbf{A}}\right)$; $\mathbf{Q} \leftarrow \hat{\mathbf{A}}^{T}\hat{\mathbf{A}}$; Update temporal: $\hat{{\textbf{\textit{c}}}}_{k} \leftarrow \max\left(0, \hat{{\textbf{\textit{c}}}}_{k} + \frac{\mathbf{P}(:,k)-\hat{\mathbf{C}}\mathbf{Q}(:,k)}{\mathbf{Q}(k,k)}\right)$; $\mathrm{corr}(k,{\textbf{\textit{x}}})\leftarrow \frac{1}{T\cdot\nu({\textbf{\textit{x}}})\cdot\mathrm{sd}({\textbf{\textit{c}}}_k)}\mathbf{U}({\textbf{\textit{x}}},:)\left((\mathbf{V} - \bar{\mathbf{V}})({\textbf{\textit{c}}}_k - \bar{{\textbf{\textit{c}}}}_k)\right)$; Update spatial support: $\Omega_k \leftarrow$ biggest connected component in $\{{\textbf{\textit{x}}}|\mathrm{corr}(k,{\textbf{\textit{x}}})\geqslant\epsilon_1\}$\
that spatially overlaps with $\{{\textbf{\textit{a}}}_k>0\}$; $\hat{{\textbf{\textit{a}}}}_k(\Omega_k^{c}) \leftarrow 0$; $\rho(k,{\textbf{\textit{x}}})\leftarrow\left(\mathrm{corr}(k,{\textbf{\textit{x}}})\geqslant\epsilon_2\right)$; Merge overlapping components $k_1,k_2$ if $\sum_{{\textbf{\textit{x}}}} \left(\rho(k_1,{\textbf{\textit{x}}}) * \rho(k_2,{\textbf{\textit{x}}})\right) / \sum_{{\textbf{\textit{x}}}}\rho(k_i,{\textbf{\textit{x}}}) \geqslant \epsilon_3$; $(\tilde{{\textbf{\textit{a}}}},\tilde{{\textbf{\textit{c}}}}) \leftarrow$ rank-1 NMF on $[\hat{{\textbf{\textit{a}}}}_{k_1},\cdots,\hat{{\textbf{\textit{a}}}}_{k_r}][\hat{{\textbf{\textit{c}}}}_{k_1},\cdots,\hat{{\textbf{\textit{c}}}}_{k_r}]$ for merged components $k_1,\cdots,k_r$; $\hat{\mathbf{A}}\leftarrow \left[\hat{\mathbf{A}}\backslash \{{\textbf{\textit{a}}}_{k_1},\cdots,{\textbf{\textit{a}}}_{k_r}\},\tilde{{\textbf{\textit{a}}}}\right], \hat{\mathbf{C}}\leftarrow \left[\hat{\mathbf{C}}^T\backslash \{{\textbf{\textit{c}}}_{k_1},\cdots,{\textbf{\textit{c}}}_{k_r}\},\tilde{{\textbf{\textit{c}}}}\right]^T;$ update number of components $K$; $\hat{\mathbf{A}},\hat{\mathbf{C}},\hat{{\textbf{\textit{b}}}}$

Results {#results .unnumbered}
=======

Denoising {#denoising .unnumbered}
---------

  -------------- -------- ------------ ----------------- ------------------ ----------- ----------------- ------------
   **Dataset**             **Method**   **Compression**      **Total**        **SNR**                     
                  Frames      FOV            Patch                           **ratio**   **runtime (s)**   **metric**
    Endoscopic     6000     256x256          16x16         Patch-wise PMD       23            220.4           2.3
                                             16x16        Patch-wise PCA\*       X              X              X
                                              NA            Standard PCA         2            595.5           1.3
    Dendritic      1000     192x192          16x16         Patch-wise PMD       52             3.2            3.7
                                             16x16         Patch-wise PCA       32             1.2            2.5
                                              NA            Standard PCA         2            18.3            1.1
   Three-photon    3650     160x240          20x20         Patch-wise PMD       94            12.4            1.8
                                             20x20         Patch-wise PCA       44             3.5            1.4
                                              NA            Standard PCA         2            187.2           1.0
    Widefield      1872     512x512          32x32         Patch-wise PMD       298           12.5            3.5
                                             32x32         Patch-wise PCA       265           10.1            3.4
                                              NA            Standard PCA        10            80.1            1.6
     Voltage       6834      80x800          40x40         Patch-wise PMD       180           30.5            2.8
                                             40x40         Patch-wise PCA       213            8.7            2.7
                                              NA            Standard PCA         8            185.1           1.0
  -------------- -------- ------------ ----------------- ------------------ ----------- ----------------- ------------

  : Summary of performance for PCA vs. PMD(TV,TF). SNR metric: average ratio of denoised vs raw SNR, with average restricted to top 10% of pixels with highest raw SNR (to avoid division by small numbers when calculating SNR ratios); an SNR metric of 1 indicates no improvement compared to raw data. Compression ratio defined in the main text. \* denotes that the patch-wise PCA method left a significant amount of visible signal in the residual for this dataset, and therefore we did not pursue further comparisons of timing or the other statistics shown here. To obtain optimistic results for the standard PCA baseline, runtimes are reported for a truncated SVD with prior knowledge of the number of components to select for each dataset (i.e., runtimes did not include any model selection steps for standard PCA). Results for patch-wise methods are reported for a single (non-overlapping) tiling of the FOV; note that total runtimes are reported (not runtimes per patch). All experiments were run using an Intel Core i7-6850K 6-core processor.[]{data-label="tab:pro_pro"}

![Illustration of the compression approach applied to microendoscopic imaging data. Top: individual frame extracted from the raw movie $\mathbf{Y}$ (left), denoised movie $\hat{\mathbf{Y}}$ (middle), and residual $\mathbf{Y} - \hat{\mathbf{Y}}$ (right). Bottom: example single-pixel traces from the movie (locations of pixels are circled in the top plots; first trace indicated by the black circle and second trace indicated by the gray circle). Note that the denoiser increases SNR significantly, and minimal signal is left behind in the residual. These results are best viewed in video form; see [microendoscopic imaging video](\VideoEndoscopeURL) for details.[]{data-label="fig:denoised_endoscope_1"}](./pmd_results/Endoscope_PMD.pdf){width="18cm" height="14cm"}

![Further analysis of microendoscopic imaging data. Top: per-pixel SNR estimated from the raw movie $\mathbf{Y}$ (left), denoised movie $\hat{\mathbf{Y}}$ (middle), and residual $\mathbf{Y} - \hat{\mathbf{Y}}$ (right). Red box indicates zoomed-in region shown in the previous figure. Bottom left panel: ratio of denoised vs. raw SNR; compression boosts SNR by roughly a factor of two here. Bottom middle and right: “correlation images" quantifying the average correlation of the temporal signals in each pixel vs. those in the nearest neighbor pixels [@Smith_2010], computed on raw and residual data, indicating that minimal signal is left behind in the residual. All results here and in the previous figure are based on background-subtracted data, for better visibility. []{data-label="fig:denoised_endoscope_2"}](./pmd_results/Endoscope_PMD.pdf){width="18cm" height="14cm"}

![Example frames and traces from Bessel dendritic imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_1\]. See [Bessel dendritic imaging demixing video](\VideoDemixDendriticURL) for details.[]{data-label="fig:denoised_dendritic_1"}](./pmd_results/Dendritic_PMD.pdf){width="18cm" height="14cm"}

![Summary quantification for denoising of Bessel dendritic imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_2\]. []{data-label="fig:denoised_dendritic_2"}](./pmd_results/Dendritic_PMD.pdf){width="18cm" height="14cm"}

![Example frames and traces from three-photon imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_1\]. See [three-photon imaging video](\VideoThreePURL) for details. []{data-label="fig:denoised_3p_1"}](./pmd_results/3P_PMD.pdf){width="18cm" height="14cm"}

![Summary quantification for denoising of three-photon imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_2\]. []{data-label="fig:denoised_3p_2"}](./pmd_results/3P_PMD.pdf){width="18cm" height="14cm"}

![Example frames and traces from widefield imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_1\]. See [widefield imaging video](\VideoWidefieldURL) for details. []{data-label="fig:denoised_widefield_1"}](./pmd_results/Widefield_PMD.pdf){width="18cm" height="14cm"}

![Summary quantification for denoising of widefield imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_2\].[]{data-label="fig:denoised_widefield_2"}](./pmd_results/Widefield_PMD.pdf){width="18cm" height="14cm"}

![Example frames and traces from voltage imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_1\]. See [voltage imaging demixing video](\VideoDemixVoltageURL) for details.[]{data-label="fig:denoised_voltage_1"}](./pmd_results/QState_PMD.pdf){width="18cm" height="14cm"}

![Summary quantification for denoising of voltage imaging data. Conventions as in Figure \[fig:denoised\_endoscope\_2\]. []{data-label="fig:denoised_voltage_2"}](./pmd_results/QState_PMD.pdf){width="18cm" height="14cm"}

We have applied the denoising and compression approach described above to a wide variety of functional imaging datasets (See Appendix for full details):

-   **Endoscopic**: one-photon microendoscopic calcium imaging in dorsal striatum of behaving mouse

-   **Dendritic**: two-photon Bessel-beam calcium imaging of dendrites in somatosensory cortex of mouse in vivo

-   **Three-photon**: three-photon calcium imaging of visual cortex of mouse in vivo

-   **Widefield**: one-photon widefield whole-cortex calcium imaging in behaving mouse

-   **Voltage**: one-photon in vitro voltage imaging under optogenetic stimulation.

The proposed methods perform well in all cases with no parameter tuning. We obtain compression ratios (defined as $nnz(\mathbf{Y}) / [nnz(\mathbf{U})+nnz(\mathbf{V})]$, where $nnz(\mathbf{A})$ counts the number of nonzero elements of the matrix $\mathbf{A}$) of 20x-200x, and SNR improvements typically in the range of about 2x but ranging up to 10x, depending on the dataset and the region of interest (we find that SNR improvements are often largest in regions of strongest activity, so SNR improvements vary significantly from pixel to pixel). See Table \[tab:pro\_pro\] and Figures \[fig:denoised\_endoscope\_1\]-\[fig:denoised\_voltage\_2\] for details.

In terms of runtime, we observed the expected scaling: the proposed method scales linearly in $T$, $d$, and the number of extracted components. In turn, the number of estimated components scales roughly proportionally to the number of neurons visible in each movie (in the datasets with single-cell resolution). Total runtimes ranged from a few seconds to a few minutes (for the “Endoscope" dataset, which had the largest number of extracted components); these runtimes are fast enough for the proposed method to be useful as a pre-processing step to be run prior to demixing.

We also performed comparisons against two simpler baselines: standard PCA run on the full dataset, and “patch-wise PCA" run on the same patches as used by PMD. For patch-wise PCA, we used the same stopping rule for choosing the rank of $\hat{\mathbf{Y}}$ as described above for PMD, but did not apply the TV or TF penalty. We find that using the same rank selection criterion for PCA applied to the full dataset performs relatively poorly: in each of the five datasets examined, this approach left significant visible signal behind in the residual. Thus, to make the comparisons as favorable as possible for standard PCA, we chose the rank manually, to retain as much visible signal as possible while keeping the rank as low as possible. Nonetheless, we found that the PMD approach outperformed standard PCA significantly on all three metrics examined here (compression ratio, SNR improvement, and runtime), largely because PCA on the full image outputs dense $\mathbf{U}$ matrices (leading to slower computation and worse noise suppression) whereas the $\mathbf{U}$ matrices output by the patch-wise approaches are highly sparse.

The patch-wise PCA approach has much stronger performance than standard PCA applied to the full data. In four out of five datasets (the “Endoscope" dataset was the exception) patch-wise PCA captured all the visible signal in the dataset and did not leave any visible signal behind in the residual. In these four datasets PMD performed comparably or significantly better than patch-wise PCA in terms of SNR improvement and compression score, but patch-wise PCA was faster. Thus there may be some room to combine these two approaches, e.g., to use PCA as a fast initial method and then PMD to provide further denoising and compression. We leave this direction for future work.

Demixing {#demixing-1 .unnumbered}
--------

### Voltage imaging data {#voltage-imaging-data .unnumbered}

  -- --
     
     
     
     
     
     
  -- --

  -- -- -- -- -- -- --
                    
                    
                    
                    
                    
                    
                    
                    
                    
                    
  -- -- -- -- -- -- --

Next we turn to the problem of demixing. We begin with an analysis of a challenging voltage imaging dataset. Voltage imaging (VI) data presents a few important challenges compared to calcium imaging (CI) data: currently-available VI data typically has much lower SNR and displays much stronger bleaching effects than CI data. The dataset we focus on here has another challenging feature: the preparation was driven with time-varying full-field optogenetic stimulation, resulting in highly correlated subthreshold activity in the visible cells, which are highly overlapping spatially. In preliminary analyses of this data we applied variants of CNMF-E [@zhou2018efficient] but did not obtain good results (data not shown), due to the strong bleaching and optogenetic stimulation-induced correlations present in this data.

Thus we pre-processed this data by applying a spline-based detrending to each pixel (see Appendix for full details). This served to attenuate the highly-correlated bleaching signals and subthreshold fluctuations in the raw data, leaving behind spiking signals (which were not perfectly correlated at the millisecond resolution of the video data here) along with uncorrelated noise as the dominant visible signals in the data. Figure \[fig:vi\_superpixels\] shows that the denoiser (followed by soft-thresholding) serves to significantly improve the separability of neural signals from noise in this data: the superpixels obtained after denoising and soft-thresholding provide excellent seeds for the constrained NMF analysis. Figures \[fig:vi\_demixing\] (and the corresponding video) and \[fig:vi\_components\] demonstrate that the full demixing pipeline achieves good performance, extracting components with high spatial and temporal SNR and leaving relatively little residual signal behind despite the limited SNR and the multiple overlapping signals visible in the original (detrended) data. Note that in the final step we project the estimated spatial components back onto the original data, recovering the (highly correlated) temporal components including strong bleaching components (panel D of Figure \[fig:vi\_components\]). Finally, we achieved a speedup in the NMF iterations here that was roughly proportional to the ratio of the rank of $\mathbf{Y}$ compared to the rank of $\mathbf{U}$.

### Bessel dendritic imaging data {#bessel-dendritic-imaging-data .unnumbered}

[ccc]{} Proposed pipeline & NMF on $\hat{\mathbf{Y}}$ & NMF on $\mathbf{Y}$\

[cccc]{} Ground truth & Proposed pipeline & NMF on $\hat{\mathbf{Y}}$ & NMF on $\mathbf{Y}$\

[ccc]{} Spatial components & Spatial components support & Temporal components\
\

The VI dataset analyzed in the preceding subsection contained a number of large visible axonal and dendritic components, but also displayed strong somatic components. For our next example we focus on a CI dataset dominated by dendritic components, where the simple Gaussian spatial filter approach introduced in [@pnevmatikakis2016simultaneous] for initializing somatic components is ineffective. (Indeed, in dendritic or axonal imaging datasets, a search for “hotspots" in the images is biased towards pixels summing activity from multiple neurons — and these “non-pure" pixels are exactly those we wish to avoid in the demixing initialization strategy proposed here.)

Figure \[realcompare\] illustrates several of the spatial components extracted by our pipeline (again, see the corresponding video for a more detailed illustration of the demixing performance); these components visually appear to be dendritic segments and match well with the signals visible in the data movie. Notably, no parameter tuning was necessary to obtain good demixing performance on both the VI and CI datasets, despite the many differences between these data types. Additionally, as a baseline comparison we applied a simple sparse NMF approach with random initialization (similar to the method described in [@pnevmatikakis2016simultaneous]) to both the denoised and raw data ($\hat{\mathbf{Y}}$ and $\mathbf{Y}$, respectively). As shown in the right columns of Figure \[realcompare\], this baseline approach extracted components that were much more mixed and noisy than the components extracted by our proposed demixing pipeline; we also found that the baseline approach was more prone to missing weaker, dimmer components than was the proposed pipeline (data not shown).

The above analyses depended on qualitative visual examinations of the obtained components and demixing video. We also generated simulated data with characteristics closely matched to the raw data, in order to more quantitatively test the demixing performance against a known (albeit simulated) ground truth. To generate simulated data $\mathbf{Y}$, we used the $\mathbf{A}$ and $\mathbf{C}$ estimated from the raw data, and further estimated the conditional distribution of the residual as a function of the denoised data $\mathbf{A} \mathbf{C}$ in the corresponding pixel $x$ and time bin $t$; then we added independent noise samples from this signal-dependent conditional distribution (but with the noise scale multiplied 2x, to make the simulation more challenging) to $\mathbf{AC}$. See the [simulated Bessel dendritic imaging video](\VideoSimulateDendriticURL) for comparison of real and simulated data. We ran the three demixing pipelines on this simulated data. Typical results of these simulations are shown in Figure \[simcompare\]: again we see that the proposed pipeline captures the ground truth components much more accurately than do the baseline methods, similar to the results shown in Figure \[realcompare\]. Quantitatively, components extracted by proposed pipeline have higher correlation with ground truth components than do those extracted by sparse NMF approaches, as shown in Figure \[corr\_sim\_comp\].

Discussion {#discussion .unnumbered}
==========

We have presented new scalable approaches for compressing, denoising, and demixing functional imaging data. The compression and denoising methods presented are generally applicable and can serve as a useful generic step in any functional video processing pipeline, following motion correction and artifact removal. The new demixing methods proposed here are particularly useful for data with many dendritic and axonal processes, where methods based on simple sparse NMF are less effective.

Related work {#related-work .unnumbered}
------------

Other work [@haeffele2014structured; @pnevmatikakis2016simultaneous; @de2018structured] has explored penalized matrix decomposition incorporating sparsity or total variation penalties in related contexts. An important strength of our proposed approach is the focus on highly scalable patch-wise computations (similar to [CaImAn](https://github.com/flatironinstitute/CaImAn)); this leads to fast computations and avoids overfitting by (implicitly) imposing strong sparsity constraints on the spatial matrix $\mathbf{U}$. We also employ a constrained optimization approach using the trend-filtering (TF) penalty, which is more flexible e.g. than the sparse convolutional temporal penalty used in [@haeffele2014structured], since the constrained TF approach doesn’t require us to fit a specific convolutional model or to estimate any Lagrange multipliers for the sparsity penalty.

There are also some interesting connections between the demixing approach proposed in [@petersen2017scalpel] and our approach to initializing NMF, which is based on the sparse projection algorithm (SPA). [@fu2015self; @gillis2018fast] discuss the relationships between SPA and group-sparse dictionary selection methods related to the approach used in [@petersen2017scalpel]; thus the methods we use to compute “pure" superpixels and the methods used in [@petersen2017scalpel] to select neural dictionary elements are closely related. However, our denoise-then-superpixelize approach to seeding the dictionary of neural temporal components is in a sense converse to the clustering approach developed in [@petersen2017scalpel] for seeding the dictionary of neural spatial components. There may be room to fruitfully combine these two approaches in the future.

Future work {#future-work .unnumbered}
-----------

Real-time online updates for $\mathbf{U}$ and $\mathbf{V}$ should be possible, which would enable the incorporation of the compression and denoising approach into [@giovannucci2017onacid] for improved online demixing of neural activity. We are also continuing to explore alternative methods for spatial and temporal denoising of $\mathbf{u}_k$ and $\mathbf{v}_k$, e.g. artificial neural network denoisers.

In the near future we plan to incorporate our code into the [CaImAn](https://github.com/flatironinstitute/CaImAn) and [CNMF-E](https://github.com/zhoupc/CNMF_E) packages for calcium imaging analysis. We hope that the proposed compression methods will help facilitate more widespread and routine public sharing of these valuable datasets and lead to more open and reproducible neuroscience.

Code {#code .unnumbered}
----

Open source code is available at [https://github.com/paninski-lab/funimag](https://github.com/paninski-lab/funimag
).

Video captions {#video-captions .unnumbered}
--------------

1.  \
    (left) Raw movie $\mathbf{Y}$; (middle) background $\mathbf{Y}_{BG}$ estimated via rank-5 PMD; (right) estimated foreground $\mathbf{Y} - \mathbf{Y}_{BG}$. Ticks along the horizontal and vertical axis (in this video and in the videos below) indicate patch borders; note that no edge artifacts are visible at these borders.

2.  \
    (left) Foreground; (middle) denoised foreground $\hat{\mathbf{Y}}$; (right) residual $\mathbf{Y} - \hat{\mathbf{Y}}$.

3.  \
    (left) Raw movie $\mathbf{Y}$; (middle) denoised movie $\hat{\mathbf{Y}}$; (right) residual $\mathbf{Y} - \hat{\mathbf{Y}}$.

4.  \
    Same format as previous video.

5.  \
    Panels from top to bottom: (1) detrended movie $\mathbf{Y}$; (2) denoised movie $\hat{\mathbf{Y}}$; (3) MAD soft-thresholded movie; (4) rank-1 NMF approximation within superpixels; (5) superpixels; (6) pure superpixels.

6.  \
    Panels from top to bottom: (1) detrended movie $\mathbf{Y}$; (2) denoised movie $\hat{\mathbf{Y}}$; (3) estimated signal $\mathbf{AC}$; (4) background $\mathbf{B}$; (5) residual $\hat{\mathbf{Y}} - \mathbf{AC} - \mathbf{B}$; (6) estimated noise $\mathbf{Y} - \hat{\mathbf{Y}}$.

7.  \
    Top: (left) motion corrected movie $\mathbf{Y}$; (middle) denoised movie $\hat{\mathbf{Y}}$; (right) estimated signal $\mathbf{AC}$; Bottom: (left) background $\mathbf{B}$, (middle) residual $\hat{\mathbf{Y}} - \mathbf{AC} - \mathbf{B}$, and (right) estimated noise $\mathbf{Y} - \hat{\mathbf{Y}}$.

8.  \
    Top: (left) Motion corrected real movie; (right) simulated movie. Bottom: (left) estimated noise from real movie; (right) simulated noise.

Acknowledgments {#acknowledgments .unnumbered}
---------------

We thank Shay Neufeld and Bernardo Sabatini for generously sharing their micro-endoscopic data with us, and Andrea Giovanucci, Eftychios Pnevmatikakis, Ziqiang Wei, Darcy Peterka, Jack Bowler, and Uygar Sumbul for helpful conversations. We also thank our colleagues in the International Brain Laboratory for motivating our efforts towards compressing functional imaging data. This work was funded by Army Research Office W911NF-12-1-0594 (MURI; EH and LP), the Simons Foundation Collaboration on the Global Brain (LP), National Institutes of Health R01EB22913 (LP), R21EY027592 (LP), 1U01NS103489-01 (NJ and LP), R01NS063226 (EH), R01NS076628 (EH), RF1MH114276 (EH), and U19NS104649-01 (EH and LP); in addition, this work was supported by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/ Interior Business Center (DoI/IBC) contract number D16PC00003 (LP). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author contributions {#author-contributions .unnumbered}
--------------------

EKB and LP conceived the project. EKB led development of the local PCA compression and denoising approach, including the 4x overcomplete approach for avoiding block artifacts. IK led development of the PMD(TF,TV) approach. DZ led development of the superpixelization and local NMF demixing approach. RZ developed a preliminary version of the PMD approach. PZ contributed to the development of the demixing approach. FG, JF, and GD contributed the voltage imaging dataset. JR, PF, TM, and AT contributed the three-photon imaging dataset. YL, RL, and NJ contributed the Bessel dendritic dataset. YM, SK, MS, and EH contributed the widefield dataset. EKB, IK, DZ, and LP wrote the paper, with input from PZ. LP supervised the project.

Appendix: dataset details {#appendix-dataset-details .unnumbered}
=========================

Microendoscopic imaging data {#microendoscopic-imaging-data .unnumbered}
----------------------------

This dataset was analyzed previously in [@zhou2018efficient]; see the “Dorsal Striatum Data" subsection of the Methods section of that paper for full experimental details. Briefly, a 1 mm gradient index of refraction (GRIN) lens was implanted into dorsal striatum of a mouse expressing AAV1-Syn-GCaMP6f; imaging was performed using a miniature one-photon microscope with an integrated 475 nm LED (Inscopix) while the mouse was freely moving in an open-field arena. Images were acquired at 30 Hz and then down sampled to 10 Hz.

Bessel dendritic imaging data {#bessel-dendritic-imaging-data-1 .unnumbered}
-----------------------------

All surgical procedures were in accordance with protocols approved by the Howard Hughes Medical Institute Janelia Research Campus Institutional Animal Care and Use Committee. C57BL/6J mice over 8 weeks old at the time of surgery were anesthetized with isoflurane anesthesia (1–2%). A craniotomy over nearly the entire left dorsal cortex (from Bregma +3 mm to Bregma -4.0 mm) was performed with the dura left intact, with the procedure described in detail previously in [@sofroniew2016large]. AAV2/9-synapsin-flex-GCaMP6s (2.5$\times 10^{13}$ GC/ml) was mixed with AAV2/1-synapsin-Cre (1.5$\times 10^{13}$ GC/ml, 1000$\times$dilution with PBS) at 1:1 to make the working viral solution for intracerebral injections. 30 nl viral solution was slowly injected into exposed cortex at 0.5 mm below dura. Injection sites were evenly spaced (at 0.7-0.9 mm separation) along two lines at 2.3 mm and 3.3 mm parallel to the midline. A custom-made glass coverslip (450 $\mu$m thick) was embedded in the craniotomy and sealed in place with dental acrylic. A titanium head bar was attached to the skull surrounding the coverslip. After recovery from surgery, the mice were habituated to head fixation. Four weeks after surgery, the head-fixed mouse was placed on a floating ball in the dark. The spontaneous neural activity as indicated by GCaMP6s fluorescence signal was recorded in the somatosensory cortex.

Volumetric imaging of dendrites was achieved by scanning an axially extended Bessel focus in [@lu201850] and [@lu2017video]. An axicon-based Bessel beam module was incorporated into a 2-photon random access mesoscope (2p-RAM) in [@lu201850]. Details of the 2p-RAM have been described previously in [@sofroniew2016large]. Briefly, the system was equipped with a 12kHz resonant scanner (24 kHz line rate) and a remote focusing unit that enabled fast axial movements of the focal plane. The system has an excitation numerical aperture (NA) of 0.6 and a collection NA of 1.0. The measured lateral full width at half maximum (FWHM) of the Gaussian focus at the center of the field of view was  0.65 $\mu$m. The lateral and axial FWHMs of the Bessel focus were 0.60 $\mu$m and 71 $\mu$m, respectively. Scanning the Bessel focus in two dimensions, therefore, probed brain volumes within a  100 $\mu$m axial range. The volumetric dendritic data presented in this paper were obtained by placing the center of the Bessel focus at 62 $\mu$m below dura to probe structures at 12 $\mu$m to 112 $\mu$m below dura (figure \[gaussian\_bessel\]). Dendrites within this volume were imaged at an effective volume rate of 3.7 Hz, with each image having 1924$\times$2104 pixels at 0.33 $\mu$m/pixel in the x-y plane. The wavelength of the excitation light was 970 nm and the post-objective excitation power was 120 mW. Images were spatially decimated and cropped for the analyses shown here.

![In vivo volumetric imaging of dendrites in the mouse brain. (a) Maximum intensity projection of a 3D volume (635 $\mu$m x 694 $\mu$m x 100 $\mu$m) of dendrites. The sampling size was 0.33 $\mu$m/pixel. Post-objective power: 24 mW. (b) Image of the same volume collected by scanning a Bessel focus with 0.60 $\mu$m lateral FWHM and 71 $\mu$m axial FWHM. The effective volume rate was 3.7 Hz. Post-objective power: 120 mW. Excitation wavelength: 970 nm. Scale bar: 100 $\mu$m. []{data-label="gaussian_bessel"}](./plots/experiment/gaussian_bessel.png){width="100.00000%"}

Three-photon imaging data {#three-photon-imaging-data .unnumbered}
-------------------------

All procedures were carried out in accordance with the ethical guidelines of the National Institutes of Health and were approved by the Institutional Animal Care and Use Committee (IACUC) of Baylor College of Medicine. Cranial window surgeries over visual cortex were performed as described previously [@reimer2014pupil]. Briefly, a 4 mm cranial window was opened under isoflurane anesthesia and sealed with a 4 mm glass coverslip and surgical glue. The dura was removed before applying the coverslip to increase optical access to the cortex. Imaging was performed in a triple-transgenic mouse (Slc17a7-Cre x Dlx5-CreER x Ai148) expressing GCaMP6f pan-neuronally throughout cortex. Three-photon imaging data was collected as described previously [@ouzounov2017vivo]. Three-photon excitation of GCaMP6 was at 1320nm, which also enabled visualization of unlabeled vasculature and white matter via THG (third harmonic generation). Power was calibrated prior to each day of scanning and carefully maintained below 1.5nJ at the focal plane. For this study, scans were collected at 680 microns and 710 microns below the cortical surface with a 540 x 540 micron field of view at 0.59 pixels/micron spatial resolution and a frame rate of 5 Hz. Imaging was performed at the border of V1 and LM during presentation of oriented noise stimuli.

Widefield imaging data {#widefield-imaging-data .unnumbered}
----------------------

See [@ma2016resting; @ma2016wide] for full details.

Voltage imaging data {#voltage-imaging-data-1 .unnumbered}
--------------------

Q-State’s proprietary Optopatch all-optical electrophysiology platform was used to record fluorescence recordings from induced pluripotent stem (iPS) cell-derived NGN2 excitatory neurons from a cohort of human subjects [@werley2017all]. Stimulation of action potentials was achieved with a blue light-activated channelrhodopsin (CheRiff). Fluorescent readout of voltage was enabled by an Archaerhodopsin variant (QuasAr). NGN2 neurons were produced at Q-State using a transcriptional programming approach. Recordings were performed with an ultra-widefield instrument with a resolution of 800x80 pixels (corresponding field of view of 2 mm$^2$) at a frame rate of 987 Hz.

The obtained data displayed breaks during stimulus resets and photobleaching. To remove these effects from the raw data, we removed frames during stimulus resets, extracted slow trends with a robust B-spline regression (with knots chosen to allow for non-differentiability at stimulus change-points and discontinuity at stimulus resets), and then a quadratic regression against frames with no stimuli to capture and then remove photobleaching effects.

[^1]: Equal contribution, arranged alphabetically; ekb2154, iak2119, dz2336@columbia.edu

[^2]: Departments of Statistics and Neuroscience, Grossman Center for the Statistics of Mind, Center for Theoretical Neuroscience, and Zuckerman Mind Brain Behavior Institute, Columbia University

[^3]: Q-State Biosciences, Inc., Cambridge, MA

[^4]: Department of Biomedical Engineering and Zuckerman Mind Brain Behavior Institute, Columbia University

[^5]: Departments of Physics and Molecular and Cell Biology, UC Berkeley

[^6]: Department of Neuroscience and Center for Neuroscience and Artificial Intelligence, Baylor College of Medicine

[^7]: One important note: many matrix factorizations are possible here to obtain a compressed representation $(\mathbf{U},\mathbf{V})$. This non-uniqueness does not pose an issue for either compression or denoising. This makes these problems inherently easier than the demixing problem, where the identifiability of $\mathbf{A}$, $\mathbf{C}$, and $\mathbf{B}$ (perhaps up to permutations of the rows and columns of $\mathbf{A}$ and $\mathbf{C}$) is critical.
