Manon Ansart
email: manon.ansart@u-bourgogne.fr

Aubert Thierry Cresson
email: thierry.cresson@etsmtl.ca

Vázquez Jacques Benjamin

Benjamin Aubert
email: baubert@eos-imaging.com

Jacques De Guise
email: jacques.deguise@etsmtl.ca

Carlos Vázquez
email: carlos.vazquez@etsmtl.ca

Statistical shape model of the spine fitting study: impact of clipping the latent representation

Keywords: Satistical shape model, principal component analysis, clipping, scoliosis, spine, vertebrae, automatic reconstruction

Statistical shape models (SSMs) represent the distribution of labeled points across a training set of shapes. The standard practice for SSMs based on principal component analysis (PCA) is to use clipping, thresholding the latent representation so that all shapes lie within 3 standard deviations of the mean. This practice precludes the representation of shapes that are not well represented by the training set, constraining the model to realistic solutions, but making it impossible to work with shapes at the edges of the statistical population. In this study, we investigate the impact of clipping in a PCA-based SSM and whether using L2 regularization is a good replacement for clipping in the context of the automatic 2D to 3D reconstruction of the spine. We first show that using L2 regularization is equivalent to using a probabilistic PCA with two error variables, accounting for the suppression of the least important principal components and for the fact that the training set cannot perfectly represent all shapes at test time. Secondly, we use two data sets of 1746 and 768 patients with adolescent idiopathic scoliosis to study the effect of regularization, for different regularization weights and with or without clipping, for removing landmark detection errors using a simulated noise or a reconstruction pipeline. In both sets of experiments, we show that regularization removes noise in a way similar to clipping without preventing the reconstruction of out-of-distribution shapes, leading to outputs closer to ground truth, demonstrating that using a regularized SSM should be preferred to clipping.

INTRODUCTION

Statistical shape models (SSM) learn the typical distribution of labeled points from a training data set. They are used to bring statistical information when detecting shapes in images, for segmentation tasks and 3D shape reconstruction from 2D images, such as in 1 for the fully automated reconstruction of the 3D scoliotic spine.

The first statistical shape model, proposed in 2, relies on an average shape and deformation modes computed using principal component analysis (PCA). This original method proposed to set limits on the latent representation, as parameters are assumed to follow a gaussian distribution and 99.7% of the data set should lie within three standard deviations (std) of the average shape. This procedure known as "clipping" does not allow the reconstruction of shapes that are poorly represented by the training set. Since then, different PCA-based SSMs have been proposed to represent vertebrae models, [START_REF] Whitmarsh | 3D reconstruction of the lumbar vertebrae from anteroposterior and lateral dual-energy X-ray absorptiometry[END_REF][START_REF] Khallaghi | Registration of a Statistical Shape Model of the Lumbar Spine to 3D Ultrasound Images[END_REF][START_REF] Kirschner | 3D Active Shape Model Segmentation with Nonlinear Shape Priors[END_REF] the whole spine [START_REF] Aubert | Toward Automated 3D Spine Reconstruction from Biplanar Radiographs Using CNN for Statistical Spine Model Fitting[END_REF][START_REF] Korez | Sparse and multiobject pose+shape modeling of the three-dimensional scoliotic spine[END_REF][START_REF] Boisvert | 3D reconstruction of the human spine from radiograph(s) using a multi-body statistical model[END_REF] or the rib cage [START_REF] Benameur | Three-Dimensional Biplanar Reconstruction of Scoliotic Rib Cage Using the Estimation of a Mixture of Probabilistic Prior Models[END_REF] for example. Several studies used variations of the PCA, using hierachical models, [START_REF] Lecron | Fast 3D Spine Reconstruction of Postoperative Patients Using a Multilevel Statistical Model[END_REF] kernel PCA, 5 a mixture of probabilistic PCAs. [START_REF] Benameur | Three-Dimensional Biplanar Reconstruction of Scoliotic Rib Cage Using the Estimation of a Mixture of Probabilistic Prior Models[END_REF] Despite the number of variations proposed around the original method, clipping remains a standard practice for 3D reconstruction, 10 but to our knowledge its impact has never been evaluated.

In 1 the authors proposed to use a regularization term which penalizes large values for coordinates in the latent space. This regularization has an effect similar to clipping, with a major difference: it attributes a higher cost to large values without making them impossible. This regularization therefore has the potential to constrain the SSM outputs to realistic shapes in a way similar to clipping while allowing the reconstruction of spine shapes poorly represented in the training set. No theoretical justification has however been proposed for this regularization, and it has never been used without clipping.

In this study, we investigate the impact of clipping in a PCA-based SSM of the spine. We use the model proposed in 1, and investigate whether the added regularization term can replace the clipping procedure, allowing to reconstruct shapes that are not well represented in the training set while correcting errors in landmark detection. Our contribution is two-fold. Firstly, we propose a theoretical justification for the PCA-based SSM with regularization proposed in 1, using a probabilistic PCA with two error variables, accounting for the suppression of the least important principal components and for the fact that the training set cannot perfectly represent all shapes at test time. Secondly, we investigate the impact of clipping this SSM in the context of the automatic 3D reconstruction of the scoliotic spine, using a simulated noise and using the reconstruction pipeline described in 1, on data sets of 768 and 1746 patients.

MATERIALS AND METHODS

Theoretical proposition

In this section, we propose a theoretical justification for the addition of a regularization term to the PCA-based SSM model proposed in 1. We first detail the solution proposed in 1 before showing that it is equivalent to maximising the posterior distribution of a probabilistic PCA model given partial data.

Existing model

Let the data set be composed of n r shapes, each one represented by landmarks contained in a row vector of size n c . The data set can be represented as a matrix X ∈ R nr×nc . When fitting the model to a new shape, only a subset of the landmarks are observed and used as targets for the model. In 1, the authors propose to minimize the distance between the target observed landmarks of a shape x o and the SSM output for coordinate i by using the objective function

m = arg min m   nc i=1 w i (x o,i -(x + Bm) i ) 2 + β 2 |m| k=1 m 2 k   (1)
B contains the representation modes, m the latent representation of the shape and x the mean shape model. This equation is equivalent to minimizing the squared distance between the observed landmarks and the PCA outputs while penalizing deviation from the mean shape using a L2 penalty, [START_REF] Hoerl | Ridge Regression: Biased Estimation for Nonorthogonal Problems[END_REF] with a regularization weight β. In addition, 1 clips the latent representation, ensuring that the output is within 3 standard deviation of the mean shape by thresholding m to a minimum of -3 and a maximum of 3. Assigning a weight of w i = 1 to all landmarks the solution to equation 1 is :

m = (B T o B o + β 2 I) -1 B T o (x o -x o ) (2)
where B o contains the rows of B corresponding to the observed landmarks. We show here that this model is equivalent to maximising the posterior distribution of a probabilistic PCA model given partial data, as proposed in 12, while accounting for two kinds of loss of information.

Maximum likelihood estimation for probabilistic PCA

Following 8, we use a probabilistic PCA model and define each shape x as a linear transformation of the latent variable m: x = x + Bm + . The latent variable m and the noise are assumed to follow gaussian distributions, such that m ∼ N (0, I d ) (I d being the identity matrix of size d × d) and ∼ N (0, σ 2 I nc ) (I nc being the identity matrix of size n c × n c ). As a result, x also follows a gaussian distribution 13 defined as x ∼ N (x, BB T + σ 2 I nc ).

The parameters B, x and σ can be chosen so as to maximize the likelihood. [START_REF] Tipping | Probabilistic Principal Component Analysis[END_REF] Noting C = BB T + σ 2 I nc and Σ the covariance matrix of X, the log-likelihood can be written as :

ln(p(X | x, B, σ 2 )) = - n r 2 ln(2π) - n r 2 ln|C| - 1 2 nr i=1 (x -x) T C -1 (x -x) = - n r 2 ln(2π) - n r 2 ln|C| -T r(C -1 Σ) (3)
where T r(A) is the trace of matrix A.

The maximum likelihood estimation of the parameters is given by

B = U d (Λ d -σ 2 I) 1 2 , σ 2 = 1 nc-d nc j=d+1 λ j , mi = (Λ d -σ 2 I) 1 2 U T d (x -x)
, where U d ∈ R nc×d contains the first d eigenvectors of the covariance matrix of X, sorted by decreasing eigenvalues, and Λ d contains the corresponding eigenvalues λ j .

Handling observed and missing landmarks

Although all coordinates are known for all shapes in the training set, at test time only a subset of landmarks detected in the images, x o are observed. The other landmarks are estimated using the statistical shape model, so as to minimize the distance on the observed landmarks. However, because the training set is not fully representative of the test set, the desired shape, using the observed landmarks, might not be perfectly represented by a linear combination of shapes in the training set. Following 12, we represent this variability by introducing a second slack variable 2 concerning only the observed landmarks:

x o = x o + B o m + o + 2
, where B o contains the rows of B corresponding to the observed entries, and o contains the corresponding elements of .

The slack variable 2 also follows a gaussian distribution 2 ∼ N (0, σ 2 2 I nc ), therefore ξ = o + 2 follows a gaussian distribution defined by ξ ∼ N (0, σ 2 ξ I nc ), with

σ ξ = σ + σ 2 . Therefore, x o | m ∼ N (x o + B o m, σ 2 ξ I q ). We note M = B T o B o + σ 2 ξ I. Following 14, the conditional probability distribution of the latent variable m knowing x o is m | x o ∼ N (M -1 B T o (x o -x o ), σ 2 M -1
). The latent vector maximising this distribution is therefore

m = (B T o B o + σ 2 ξ I) -1 B T o (x o -x o ) (4)
This solution is equivalent to the one proposed in 1 and described in equation 2, with β = σ 2 ξ .

Interpretation

Two slack variables are introduced in this model. The first one, , accounts for the loss of information due to the use of the d first eigenvectors only. The second one, 2 , accounts for the representation of shapes that do not belong to the training set. The combination of these slack variables explain the addition of the regularization term in equation 1.

Data set

Two data sets are used to perform two sets of experiments, described in section 2.3. Firstly, the AIS (adolescent idiopatic scoliosis) data set is composed of 768 AIS patients from the Sainte-Justine mother and child university hospital center, affiliated with the University of Montréal in Canada. Secondly, the mixed data set combines adults and AIS patients in a training set of 1536 patients and an evaluation set of 210 patients. All data sets were retrospectively collected after ethical approvals from ethical committees of the Montreal University Hospital Center (CHUM) and the École de technologie supérieure ( ÉTS), Montréal, Canada. They contain two low dose x-ray images (in sagittal and frontal views) acquired with the EOS system. 15 3D spine reconstruction was performed for each pair of images using the semi-automated method described in 16. For the evaluation set, reconstruction was performed by an expert before being corrected twice. The 3D spine is represented using a parametric spine model, 1 resulting in 1124 parameters including the 3D position of the pedicles, endplates and centers of each vertebra. As in 1 the origin is the middle of the two femoral head, the z axis is the axis joining the origin to the topmost vertebral center, the y axis is the axis joining the origin to the left femoral head and the x axis is perpendicular to the z and y axis and and pointing towards the front of the subject.

Experiments

Two sets of experiments were performed. In a first set of experiments, we evaluate the impact of clipping in a controlled simulated setting, by using the SSM to correct a simulated noise on vertebral centers. In a second set of experiments we evaluate the impact of clipping the SSM in a real pipeline, which aims at reconstructing a 3D model of the spine from 2D images. Both sets of experiments are described in more details bellow.

For all experiments, the SSM was evaluated by studying the 3D distance between estimated landmarks and true values for vertebral centers, pedicles and endplate centers. When PCA is performed, 99.999% of explained variance is kept.

Using simulated errors

The first set of experiments aims at evaluating the capacity of the SSM to correct simulated errors on vertebral centers by regularizing the global spine shape. This set of experiments was performed on the AIS data set by adding a simulated noise on vertebral centers for all vertebrae.

The noise added to the vertebral centers used as predictors was simulated using a gaussian mixture of two components, representing a standard gaussian noise (standard deviation of 3 mm across all vertebral levels) applied to all coordinates, and an aberrant noise (20 mm with a probability of p = 2 46 ) on the x and y components only. As sagittal and frontal views are combined, detection error in z rarely happens, 1 so no aberrant noise was added on the z component. The resulting probability function for the simulated noise of the x and y components is:

P (e) = (1 -p)N (0, σ) + p 2 N (20, σ) + p 2 N (-20, σ) (5)
As different errors can be expected for each coordinate and vertebra levels, different σ are used. The standard deviation were computed so that they are proportional to the global root mean square error reported in 1, with a standard deviation of 3 mm across all axes and vertebral levels. The combination of the standard and aberrant errors leads to an average 3D residual of 4.9 mm (standard deviation of 5 mm).

The noisy vertebral centers were the only predictors in the SSM, which was used to estimate the other spine parameters (providing pedicle and endplate landmarks) and denoise the vertebral centers. Different SSMs were considered, by varying the value of the regularization weight β between β = 10 -5 and β = 1.0, with and without clipping. 10-fold cross-validation was used to get an estimation on the whole data set.

Using the reconstruction pipeline

The second set of experiments was performed in actual conditions, on the mixed data set of adolescent and adult patients, using the automatic reconstruction method of the 3D spine from a pair of x-ray images proposed in 1.

This method combines the detection of vertebral landmarks (vertebral centers in a first stage, pedicles and endplate centers in a second stage), performed iteratively for each vertebra, with a correction by the spine SSM, performed after each vertebra detection. The whole spine is initialized and updated after each vertebra detection using the SSM. The landmark detection is performed using a convolutional neural network (CNN), which computes the optimal landmark displacements given the previous positions. In this method, the purpose of the SSM is two-fold: firstly, the SSM gives the position of landmarks which have not yet been detected and the value of spine parameters which are not computed from the image; secondly, the SSM corrects the position of landmarks that have already been detected, taking into account the global spine shape so as to output realistic spines. Four sets of parameters were considered: SSM with clipping, using β = 0.01 (original method 1 ) and β = 0.2; and SSM without clipping, using β = 0.01 and β = 0.1. The full pipeline was trained on the training set of 1536 subjects, and the evaluation of 210 subjects was used to evaluate the SSM. 1: Simulated noise: 3D distance (in mm) between output and true spine on several vertebra landmarks, across all subjects and vertebrae, for different regularization weights β, with and without clipping. Measures are given as mean (standard deviation).

RESULTS

Using simulated error

By varying the regularization weight β in the SSM model, with and without clipping, we can better understand the effect of regularization and clipping as well as the interaction between them. Figure 1 shows the impact of these parameters on the distance between the outputs and the noisy spines, used as predictors for the SSM, as well as the distance between the outputs and the true spines, which should ideally be minimized.

First looking at the results without clipping (in blue), we can see that a regularization weight close to 0 leads to no noise correction, resulting in an output equal to the noisy input and a large output to true spine distance. As β increases, the outputs get closer to the true spines and farther from the noisy ones, showing the expected noise correction. The distance to true spine reaches a minimum for a regularization weight of 0.05. Increasing β above this value leads to a slight increase in the distance to true spine, showing that the regularization results in a loss of information as well a noise correction. The improvement on noise correction is however still greater than the loss of information, until larger values of β (around 0.3). Secondly, looking at the results with clipping (in red), we can see that a regularization weight close to 0 leads to a large distance both to true and noisy spines, showing that shape information is not well represented. As β increases up to 0.01, the output gets closer to the true spine, but the distance to noisy spine also decreases. Clipping prevents the representation of global shapes in the test set by thresholding some latent coordinates. As β increases, high values are penalized, encouraging spreading across all modes, limiting the impact of thresholding the latent representation and hence of clipping. For values higher than 0.01, increasing β leads to effects similar to those observed without clipping, first removing noise, then removing shape information as well when β is too large. The loss of shape information however happens at 0.05, which is larger than without clipping, probably because increasing regularization also mitigates the negative impacts of clipping more strongly. The distance to true input obtained with clipping is however never smaller than the distance to true input obtained without clipping.

Table 1 shows the distance between the SSM output and true spine for vertebral centers, pedicles and endplate centers, with and without clipping, for different values of β. Together with figure 1, these results show that clipping has a negative impact, regardless of the choice of regularization weight, as the model without clipping removes the simulated noise (albeit not completely), while loosing less information, leading to outputs closer to the true spines for all β values. The choice of β highly depends on the noise magnitude, which was arbitrarily set here. Values from 0.05 to 0.2 seem reasonable choices, both with and without clipping, although validating this parameter properly would require access to the true noisy inputs, corresponding to the vertebral center detection in the images by the CNN. This information is not available, but experiments using the reconstruction pipeline allow more insight.

Using the reconstruction pipeline

Table 2: Reconstruction pipeline: 3D distance (in mm) between output and true spine on several vertebra landmarks, across all subjects and vertebrae, for different regularization weights β, with and without clipping. Measures are given as mean (standard deviation).

Table 2 shows the 3D residuals obtained on vertebra landmarks, with and without clipping, for different β values. For β = 0.01, the table shows that the results are better on average without clipping (1.9 mm for vertebral centers, standard deviation of 3.2 mm) than with clipping (2.1 mm, standard deviation of 3.7 mm), for all studied landmarks.

Although the average improvement is small, study of 3D residuals for each subject shows the impact of clipping more clearly. Figure 2 shows the average 3D residual across vertebral centers for each subject, with and without clipping, for β = 0.01. Most points are on or close to the diagonal line, representing subjects for whom the 3D residual is similar with and without clipping. Subjects lying exactly on this diagonal are subjects for which clipping has no effect, because no coordinates are saturated in the latent representation. Subjects at the bottom of the figure, bellow the diagonal, are subjects for whom the SSM without clipping yields better results than with clipping. These are subjects with unusual spine shapes, for which some coordinates are saturated. Figure 3 shows an example of such subjects. Figure 2 also shows about 12 subjects who have high residuals with clipping and for whom the residual stays high when removing clipping. Further examination of these subjects shows that the computed global spine curve is right, but the vertebrae are not well placed along the spinal curve. In particular, some vertebrae seem to not be detected at all, leading to a shift in a series of vertebrae which are identified as the next one. An example is provided in Figure 4. Subjects with at least one vertebral center with a z absolute error of more than 8.6 mm (minimum difference in the z coordinate between 2 adjacent vertebrae in the data set) are colored in orange in Figure 2. This figure shows that all the subjects with a high error that is not corrected when removing clipping belong to this category. Notably, the only subject with an error smaller with clipping than without belongs to this category. This issue is probably due to an error in the CNN detection which is not modeled by our simulated noise. Given that all the subjects that still have a high error without clipping belong to this category, this issue seems of particular importance and requires further study. Figure 3: Example of a subject for which the reconstruction is worse using clipping (left, average 3D residual across all vertebral centers of 4.9 mm) than without clipping (right, average 3D residual across all vertebral centers of 2.9 mm) for β = 0.01. For the true spine vertebral centers are displayed using squares linked by dashed lines. For the output vertebral centers are displayed using dots linked by plain lines.

Regarding the regularization weight, Table 2 shows that increasing β leads to worse results, with clipping (3D residual of 2.1 mm (std=3.7) for β = 0.01 and of 2.8 mm (std=3.7) for β = 0.2) and without (3D residual of 1.9 mm (std=3.2) for β = 0.01 and of 2.3 mm (std=3.1) for β = 0.1).

DISCUSSION AND CONCLUSION

We studied the impact of clipping in a SSM fit process from constraints that were simulated or detected in actual X-ray images. Both sets of experiments show that clipping limits model performance, by not allowing spine shapes that are not well represented in the training set, and therefore making the model depend heavily on the choice of training set. The regularizing effect of clipping, which corrects landmark detection errors, is not necessary when using a regularization term in the objective function used to fit the SSM on partial data. This regularization term penalizes deviation the mean shape, attributing a higher cost to shapes farther from the average shape. We showed that it is equivalent to using a slack variable accounting for the non-representativity of the training set, which attributes a higher cost to shapes not represented in the training set without making the reconstruction of these shapes impossible. We showed that this regularization term lessens the negative impact of clipping, but it also makes clipping unnecessary by correcting landmark detection errors. Because using the regularized SSM model also allows the reconstruction of out-of-distribution spine shapes, resulting in better reconstruction performance, it should be used instead of and without clipping. It is to be noted that these conclusions are only valid for fully automatic reconstruction. When the user is allowed to edit the model by manually moving landmarks, clipping restricts manual modifications to valid shapes. If clipping is removed, another method for limiting user modifications has to be used.

A downside of the addition of the regularization term is the introduction of the new hyper-parameter, β, which depends on the input errors that the SSM should remove, and should therefore be validated for each pipeline and data set. Although regarding clipping the same conclusions can be drawn from the experiments performed on simulated errors and on the reconstruction pipeline, the same cannot be said for the choice of the regularization weight. This discrepancy can have several explanations. It can first be a sign of the sensitivity of the hyperparameter, which could highly depend on the context and data. Secondly, the simulated error may not fully represent the detection errors made by the CNN. Improving error simulation would require more information on these detection errors. Lastly, besides this difference between real and simulated errors, the experimental conditions are very different: in the experiments with simulated errors the SSM is only used once, using all noisy vertebral centers as input. By contrast, in the reconstruction pipeline, the SSM is applied iteratively, after detection of each new vertebral center, and the process is repeated for pedicles and endplate centers after all vertebral centers are detected. The use of these landmarks in addition to vertebral centers during the second stage may also explain the lower error in the reconstruction pipeline compared to the experiments on simulated detection errors.

Figure 1 :

 1 Figure 1: Effect of the regularization weight β, with and without clipping, on the AIS data set with simulated errors. 3D distance (mm) is averaged over all vertebral centers over 768 patients after cross-validation.

Figure 2 :

 2 Figure 2: Mean 3D residual (mm) across vertebral centers for each of the 210 subjects in the mixed evaluation set, with and without clipping, (β = 0.01). The diagonal line represents equality with and without clipping, points bellow the line represent subjects for which the 3D residual is lower without clipping.

Figure 4 :

 4 Figure 4: Example of a subject with a shift in the vertebra labelling, with clipping (left) and without clipping (right) for β = 0.01 (average 3D residual across all vertebral centers of 11.7 mm for both). For the true spine vertebral centers are displayed using squares linked by dashed lines. For the output vertebral centers are displayed using dots linked by plain lines. The shift is highlighted using arrows.

Proc. of SPIE Vol. 12032 1203212-1 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-2 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-3 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-4 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-5 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-6 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-7 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-8 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 12032 1203212-9 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 13 Jun 2022 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

ACKNOWLEDGMENTS

We thank Dr Christine Bakhous and Dr Manuela Pacheco from EOS Imaging Inc. for performing some of the experiments and providing their expertise regarding the existing models. This work is supported in part by the Collaborative Research and Development program from NSERC and in part by EOS Imaging Inc.


João F C Mota
email: jmota@isr.ist.utl.pt

João M F Xavier
email: jxavier@isr.ist.utl.pt

Pedro M Q Aguiar
email: aguiar@isr.ist.utl.pt

Markus Püschel
email: pueschel@ece.cmu.edu

DISTRIBUTED ALGORITHMS FOR BASIS PURSUIT

Keywords: Basis pursuit, distributed algorithms, dual methods, subgradient method, method of multipliers, sensor networks

The Basis Pursuit (BP) problem consists in finding a least 1 norm solution of the underdetermined linear system Ax = b. It arises in many areas of electrical engineering and applied mathematics. Applications include signal compression and modeling, estimation, fitting, and compressed sensing. In this paper, we explore methods for solving the BP in a distributed environment, i.e., when the computational resources and the matrix A are distributed over several interconnected nodes. Special instances of this distributed framework include sensor networks and distributed memory and/or processor platforms. We consider two distribution paradigms: either the columns or the rows of A are distributed across the nodes. The several algorithms that we present impose distinct requirements on the degree of connectivity of the network and the per-node computational complexity.

INTRODUCTION

Representing a given signal as a sparse linear combination of elementary signals, drawn from an over-complete dictionary, is a problem of great theoretical and practical importance [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]. Formally, given a vector b ∈ R m , which may represent a signal of interest, and a matrix A ∈ R m×n with m < n, whose columns form an overcomplete dictionary of elementary functions, the task of expressing b as a linear combination of a small subset of columns of A can be approximated by the convex problem minimize

x 1 subject to Ax = b , (BP) where x ∈ R n is the variable to optimize and x 1 = |x1| + . . . + |xn|. Hereafter, we assume that A has full rank, in order to secure the feasibility of (BP) for any given b. Problem (BP) is known as the basis pursuit problem [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. Problem statement. The goal in this paper is to devise methods that solve (BP) when the matrix A is not available in a single node, but is rather spread among several nodes, each one having some computational power. We consider that either the columns or the rows of A are distributed across the nodes. This leads us to two problem statements: 1) (resp. 2)) Given P nodes, each one storing a subset of columns (resp. rows) of A, find network topologies along with distributed algorithms for solving (BP). In problem 1) we say that the partition of A is horizontal, while in problem 2) it is vertical.

Related work. To the best of our knowledge, there is no previous work concerning solving (BP) in a distributed environment, as studied herein. It is well known that (BP) can be reformulated as a linear program (LP) [START_REF] Boyd | Convex Optimization[END_REF] and the work in [START_REF] Yarmish | A distributed implementation of the simplex method[END_REF] proposed a generalization of the LP-simplex method for a distributed scenario, assuming a horizontal partition of A. However, the resulting algorithm requires a fully connected topology in the sense that, at each iteration, every node must be able to communicate with all the remaining ones. In this paper, we present methods with less demanding graph topologies (moreover, they are not based on LP reformulations of (BP)).

HORIZONTAL PARTITION

In this section, we consider that the columns of A are distributed among P nodes. The pth node stores the pth submatrix1 in the hor-

izontal partition A = [ A1 • • • Ap • • • AP ], where Ap ∈ R m×np , n1 + • • • + nP = n.
The P nodes have to be connected (we seek sparse network topologies) and cooperate in order to solve (BP). The methods we propose are based on duality [START_REF] Bertsekas | Nonlinear Programming[END_REF]. It is well known that, under the appropriate conditions (e.g., separability of cost function and constraints), dual programs can be tackled by distributed methods. We propose two equivalent reformulations of (BP) whose associated dual programs are amenable to distributed solution methods. These are explained in the two remaining subsections.

Bounded BP

Let R > 0 be a strict upper bound on the ∞ norm of any solution of (BP). Then, (BP) is equivalent to the bounded basis pursuit (BBP) problem minimize

x 1 subject to Ax = b x ∞ ≤ R . (BBP)
Before we proceed, we note that such an R is easily found: given any x such that Ax = b, pick R > n x ∞. Indeed, if x denotes a solution of (BP) then

x ∞ ≤ x 1 ≤ x 1 ≤ n x ∞ < R.
Lp(λ) = inf{ xp 1 -λ Apxp : xp ∞ ≤ R}.
We used the decomposition x = (x1, . . . , xp, . . . , xP ) where xp ∈ R np . It can be shown that the dual of (BP) is not unconstrained. This is the motivation for inserting the redundant constraint x ∞ ≤ R: it makes the corresponding dual unconstrained, and therefore directly applicable to a distributed setting (since the cost function is already separable). Solving the dual. We propose to solve (DBBP) by a standard subgradient method2 [START_REF] Bertsekas | Nonlinear Programming[END_REF], i.e.,

λ k+1 = λ k + α k g(λ k ), (1)
where k denotes the iteration number, α k > 0 is a stepsize and g(λ k ) is a supergradient3 of the concave dual function L at the point λ k . This algorithm is globally convergent, under an appropriate choice for the stepsizes, see [START_REF] Bertsekas | Nonlinear Programming[END_REF]. Using subgradient calculus [START_REF] Boyd | Subgradients, notes for ee364b[END_REF][START_REF] Bertsekas | Nonlinear Programming[END_REF] it can be shown that, for any given λ, the supergradient g(λ) can be chosen as g(λ) = b -P p=1 Apx p (λ) where the ith component of x p (λ) is φ(ri); here, ri denotes the inner-product between λ and the ith column of Ap and

φ(t) =    -R, if t ≤ -1 0, if -1 < t < 1 R, if t ≥ 1 .
Network topology. We propose to implement the subgradient method in a network where the P nodes are connected to a central node, see figure 1(a). The central node implements the recursion (1) as follows: (i) λ k is broadcasted to the P nodes, (ii) the pth node computes x p (λ) and sends Apx p (λ) to the central node -this is done in parallel for all P nodes -, (iii) the central node adds the received data to obtain g(λ k ) and moves to λ k+1 . It is assumed that b is available at the central node. Note that the computational complexity associated with each of the P nodes is very low.

Solving the primal. Let λ be a solution of (DBBP). Through the KKT optimality conditions of (BBP), we can infer valuable information about the solutions of (BBP) (hence of (BP) since they are equivalent). More precisely, let the support of a vector be the indexes of its nonzero entries, i.e., supp x = {i : xi = 0}. Then, it can be seen through the KKT system that, for any solution x of (BBP), we have supp x ⊂ I := {i :

|a i λ | ≥ 1} ( 2
)
where ai is the ith column of A (actually, |a i λ | > 1 cannot occur due to our assumption on R, but we omit this technical discussion).

In other words, once the central node knows λ it has access to an overestimate of the support of a solution of (BBP). In particular, for i ∈ I , we have x i = 0. Thus, at this point, the central node could solve a reduced basis pursuit problem, i.e., a smaller (BP) where the columns of A outside I are removed. Since the solutions of (BP) are expected to be highly sparse, this could represent a significant reduction in the problem dimensionality.

The aforementioned strategy relies on the knowledge of λ . In practice, the subgradient method is stopped after a finite number of iterations at the central node and only an estimate λ of a solution λ is available there. We adapt the previous scheme as follows. Let Î = {i : |a i λ| ≥ 1 -ξ} where ξ > 0 denotes a user-defined tolerance. Let A Î be the columns of A indexed by Î (the notation x Î follows a similar definition). Two cases are now possible, based on the feasibility of the linear system A Î x Î = b. Case 1: system is feasible -we let x Î be a solution of a reduced (BP) with A replaced by A Î . Case 2: system is unfeasible -we let x Î be a solution of the least-squares problem

minimize A Î x Î -b 2 x Î .
In either cases, we set x i = 0 for i ∈ Î. Case 2 is usually a consequence of a premature stop of the subgradient method.

Recall that A is stored in the P nodes, not in the central node. However, A Î is needed there to address either case 1 or 2. This information can be obtained as follows: after stopping the subgradient method, the central node sends λ to the P nodes which, in an obvious and parallel way, can compute their corresponding fragments of Î and A Î and report them back to the central node. Overall algorithm. We gather our results in the form of an overall algorithm. Clearly, we have two phases. In the first phase, the central node runs the subgradient method (1) to obtain λ. This is done by exchanging several messages with the P nodes, which work in parallel. In the second phase, λ is sent to the P nodes which report back their parts of Î and A Î . Then, the central nodes goes into case 1 or 2, and builds x .

• • • • • • central P 1 2 (a) Algorithm 1 • • • • • • φ 1 φ2 φ P P 1 2 x b (b) Algorithm 2

. , P do

Computation at node p:

• Compute x p (λ k ) • Return Apx p (λ k ) to the central node (c) Compute g(λ k ) = b -P p=1 Apx p (λ k ) (d) λ k+1 = λ k + α k g(λ k ) 2. λ = λ k+1

Augmented BP

In this subsection, we need an assumption whose utility will became clear soon. We assume that in the horizontal partition

A = [ A1 • • • Ap • • • AP ] there holds rank Ap = np, (3)
that is, each Ap has full column rank. Note that this automatically imposes that np ≤ m for all p.

The starting point for our new approach is based on a reformulation of (BP). Namely, for any given ρ > 0, we note that (BP) is equivalent to minimize

x 1 + ρ 2 Ax -b 2 subject to Ax = b (ABP)
since we just augmented the cost function with a term which evaluates to zero (because the constraints Ax = b are in force).

Dual problem. Although (BP) and (ABP) are equivalent, their duals are distinct. In particular, the dual of (ABP) is unconstrained in the sense that it corresponds to maximize L(λ) λ

(DABP) and L(λ) is finite everywhere, i.e., L(λ) > -∞ for all λ ∈ R m . Here, L(λ) = inf {L(x, λ) : x ∈ R n } and L(x, λ) = x 1 -λ (Ax -b) + ρ 2 Ax -b 2 . ( 4
)
Solving the dual. We propose to solve (DABP) through the iterations

λ k+1 = λ k + ρg(λ k ) (5)
where g(λ k ) is a supergradient of the concave dual function L(λ) at the point λ k : it can be computed as

g(λ k ) = b -Ax (λ k ) where x (λ) is a solution of minimize L(x, λ) x . (6)
The iterative scheme ( 5) is nothing else but the method of multipliers [START_REF] Bertsekas | Nonlinear Programming[END_REF] applied to (BP). Also, L(x, λ) in ( 4) is known as the augmented Lagrangian of (BP). The method of multipliers enjoys an important property: any cluster point of the sequence {x (λ k ) : k = 1, 2, . . .} solves (BP), see [START_REF] Bertsekas | Nonlinear Programming[END_REF].

Network topology. We propose to implement the method of multipliers (5) in a ring-shaped network, see figure 1(b). Note the absence of a central node. We arrive at this network configuration by reasoning as follows.

Let λ be fixed and focus on solving [START_REF] Boyd | Subgradients, notes for ee364b[END_REF] with respect to x = (x1, . . . , xp, . . . , xP ), xp ∈ R np . As we have P nodes available, we would like to distribute this task over them. Note that L(•, λ) is not separable with respect to the xp's, a well known drawback of the augmented Lagrangian approach. That is, a straightforward decoupling into P independent subproblems is not possible. We tackle the minimization of L(•, λ) by the nonlinear Gauss-Seidel method [5, page 267] (another viable method would be the diagonal quadratic approximation method [START_REF] Zakarian | Nonlinear Jacobi and -relaxation methods for parallel network optimization[END_REF], whose convergence properties can be established for this scenario, see [START_REF] Mota | Distributed algorithms for sparse approximation[END_REF]). This consists in minimizing L(•, λ) through sweeps. Each sweep updates the xp's by minimizing L(•, λ) with respect to each xp while keeping fixed xq, q = p. A more precise description follows: Nonlinear Gauss-Seidel method. Input: λ and initial point x 0 = (x 0 1 , . . . , x 0 p , . . . , x 0 P ). Output: x (λ) solving (6) 1. Initialize t = 0 2. Iterate on t until some stopping criterion is met:

For all p = 1, . . . , P find

x t+1 p = arg min xp L(x t+1 1 , . . . , x t+1 p-1 , xp, x t p+1 , . . . , x t P , λ) (7)
3. x (λ) = (x t 1 , . . . , x t p , . . . , x t P ) It is because assumption (3) holds that each minimization in ( 7) is well-posed, i.e., there is an unique minimizer. Regarding the convergence of the nonlinear Gauss-Seidel method in our context, the optimality of any cluster point of the sequence {(x t 1 , . . . , x t p , . . . , x t P ) : t = 1, 2, . . .} is guaranteed despite the non-smoothness of the cost function in [START_REF] Zakarian | Nonlinear Jacobi and -relaxation methods for parallel network optimization[END_REF], see [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF].

Consider now the minimization problem in [START_REF] Zakarian | Nonlinear Jacobi and -relaxation methods for parallel network optimization[END_REF]. It is natural to assign this task to node p. In the following, we show that it can be so. We start by noting that (7) boils down to

minimize 1 2 Apxp -γ t p 2 + 1 ρ xp 1 xp , (8)
where

γ t p := b - i<p Aix t+1 i - i>p Aix t i + 1 ρ λ.
Since node p has access to Ap, problem (8) can be tackled there if γ t p is available. Now, note the identity

γ t p = γ t p-1 -Ap-1x t+1 p-1 φp-1 +Apx t p . (9)
Here, φp-1 is available at the (p-1)th node while Apx t p is available at the pth node (from the previous iteration). Thus, at the pth node, it suffices to receive φp-1 from the (p -1)th node to build γ t p and tackle [START_REF] Mota | Distributed algorithms for sparse approximation[END_REF], e.g., through the method in [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF]. This flow of information can be implemented in a circular network like the one in figure 1(b).

So far, only problem (6) has been discussed. However, the whole method of multipliers (5) can be implemented in the proposed network by letting node 1 do the bookkeeping of λ k . Indeed, for a given λ k , suppose T sweeps around the network solve (6) for λ = λ k (one sweep corresponds to one t-iteration in the nonlinear Gauss-Seidel method). After the T sweeps, we have x t p x t+1 p (since we assume that the nonlinear Gauss-Seidel method has converged), and thus node 1 can calculate

g(λ k ) = γ T 1 -A1x T 1 - 1 ρ λ k .
This enables node 1 to update λ according to [START_REF] Bertsekas | Nonlinear Programming[END_REF]. Further, this node can trigger a new set of sweeps to minimize (6) for λ = λ k+1 because it can incorporate this knowledge into γ1 by making

γ 0 1 = γ T 1 - 1 ρ (λ k -λ k+1 ).
Solving the primal. After the dual problem (DABP) is solved, it suffices to do a cycle through the network to collect x(λ ) which corresponds to the last solutions of (7) stored in each node.

Overall algorithm. Leaving out the details of the (already seen) updates of λ, we summarize our results in the following algorithm. Algorithm 2. (Procedure for nodes p = 1, . . . , P ) Input: b at node 1 (and matrices Ap stored at each node). Output: x at node 1

1. Receive φp-1 from node p -1 2. if p = 1, t ← t + 1, unless some stopping criterion is met (in this case t ← 0) 3. Set γ t p = φp-1 + Apx t p 4. Find x t+1 p = arg minx p 1 2 Apxp -γ t p 2 + 1 ρ xp 1 5. Send φp = γ t p -Apx t+1 p to node p + 1
In practice, it is better to also update ρ in each iteration of the method of multipliers as ρ k+1 = cρ k where c ≥ 1, see [START_REF] Bertsekas | Nonlinear Programming[END_REF]. This can be easily incorporated in our algorithm (details are skipped). We used this latter version in our simulations.

Finally, notice that, due to the fact that in this algorithm only one node "works" at each iteration, it is possible to solve P simultaneous problems (with the same matrix A), using the architecture of figure 1(b).

VERTICAL PARTITION

In this section, we consider that the rows of A are distributed among the P nodes. The pth node stores the pth submatrix in the vertical partition

A =         A1 . . . Ap . . . AP        
where Ap ∈ R mp×n , m1 + • • • + mP = m. This scenario arises when a network of interconnected sensors, where each sensor only measures bp := Apx ∈ R mp , must reconstruct the sparse vector x without the nodes exchanging information about their sensing devices (rows of matrix A). Network topology. Contrary to the approaches in section 2, we will not impose any particular network structure. In other words, we let the network topology be given to us. We just assume that the given network of P nodes is connected. Let E = {(i, j)} denote Fig. 2. Example of a network in which algorithm 3 can run. The set of edges is E = {(1, 2), (2, 3), [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Yarmish | A distributed implementation of the simplex method[END_REF], [START_REF] Yarmish | A distributed implementation of the simplex method[END_REF][START_REF] Bertsekas | Nonlinear Programming[END_REF], [START_REF] Yarmish | A distributed implementation of the simplex method[END_REF][START_REF] Boyd | Subgradients, notes for ee364b[END_REF], [START_REF] Bertsekas | Nonlinear Programming[END_REF][START_REF] Boyd | Subgradients, notes for ee364b[END_REF], (6, 7)}.

the set of edges of the given network. We adopt the convention that i < j for each edge (i, j) ∈ E. See figure 2 for an example. Cloned BP. To design a distributed algorithm for this environment, we start by reformulating (BP) by cloning the optimization variable. More precisely, problem (BP) is equivalent to minimize P p=1

1 P xp 1 + ρ 2 (i,j)∈E xi -xj 2 subject to Apxp = bp, p = 1, . . . , P xi = xj, (i, j) ∈ E , (
+ (i,j)∈E λ (i,j) (xi -xj) + ρ 2 xi -xj 2 . ( 10
)
Solving the dual. The dual problem (DCBP) can be solved by the method of multipliers, see a related discussion in subsection 2.2. In the current setup, this corresponds to the iterations

λ k+1 (i,j) = λ k (i,j) + ρ x i {λ k (i,j) } -x j {λ k (i,j) } ( 11
)
where {x p ({λ (i,j) })} minimize [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF] with respect to {xp} subject to Apxp = bp, for λ (i,j) fixed at λ k (i,j) . We employ the nonlinear Gauss-Seidel method to find x {λ (i,j) } = x 1 ({λ (i,j) }), . . . , x P ({λ (i,j) }) .

After straightforward calculations, it can be shown that updating the variable xp ∈ R n inside a nonlinear Gauss-Seidel sweep amounts to solving at the pth node the program

minimize 1 P xp 1 + j∈Np (sign(j -p) λ (p,j) -ρxj) xp + Npρ k 2 xp 2 subject to Apxp = bp (12)
with respect to xp ∈ R n , where Np is the set of neighbors of node p and Np its cardinality. According to the canonical script of the nonlinear Gauss-Seidel method [START_REF] Bertsekas | Nonlinear Programming[END_REF], the minimizations are to be performed in the sequential order p = 1, 2, . . . , P ; thus, xj in (12) stands for x t+1 j if j < p and for x t j if j > p. Actually, the order of the minimizations can be arbitrary as long as there is an integer M such that each xp is updated at least once in every group of M contiguous iterations [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]. Hence, we don't need to worry about the sequence of our minimizations: we can let that sequence be determined by the topology of the graph. This also means that the algorithm is robust to instantaneous or even permanent link failures, provided the resulting graph remains connected.

Looking at (12), we see that the pth node only needs the variables x t+1 j and x t j that belong to its neighbors. Further, this knowledge is also sufficient to update, after each iteration k of the method of multipliers (which correspond to several sweeps of the nonlinear Gauss-Seidel method) the relevant subset of the dual variables {λ (i,j) }. In conclusion, this algorithm is totally distributed and can run on any connected graph.

Solving the primal. After the dual problem (DCBP) is solved, a solution to (BP) is available at any node: it corresponds to x {λ (i,j) } .

Overall algorithm. We summarize our findings. 2. else, in the idle times of the links, exchange convergence information with the neighbors; and, whenever the nonlinear Gauss-Seidel sweeps stop, update the dual variables λ (p,j) , for all j ∈ Np.

EXPERIMENTAL RESULTS

To test the behavior of the distributed algorithms in practice, we simulated them in an ordinary (single) processor, using Matlab. The results presented in table 1 validate the fact that all the algorithms converge, in general, to a point close to an optimal solution, as predicted by the theory. However, as each algorithm has more than one parameter to tune, we cannot draw strict conclusions with respect to the running times: we would need real distributed simulations, varying exhaustively the parameters. Since the paradigms for the horizontal and vertical partitions are different, we cannot directly compare the results between the respective algorithms.

The first two rows of table 1 show the average relative error in the variable x and in the cost function f (x), respectively. The correct solutions, x and f (x ), were considered to be the output of a Matlab linear program solver that solved (BP) in its linear program format. The numbers of the rows that contain time measurements are the average of the total time that the algorithms took to solve each problem, the respective time spent in operations that could be done in parallel (note that these simulations were performed on a single processor), and the estimated time for the execution of the algorithms if we really had P nodes available. Note that the timings are for illustration only as they do not take into account communication time or other effects of a real distributed environment. Horizontal partition. For the horizontal partition, each node (P = 10) stored 25 columns of a matrix A ∈ R 50×250 , where each entry of the matrix was drawn independently from a Gaussian distribution. Obviously, the simulated environments were those of figure 1. We first ran algorithm 2 and, after that, we ran algorithm 1 stopping it when the dual function achieved 99% of the value returned by algorithm 2. While the errors in x and f (x) for the algorithm 2 are representative of what happened in all experiments (never achieving exactly the correct solution), the same errors for the algorithm 1 are not. In fact, in more than half of the experiments, algorithm 1 returned the exact optimal solution. Concerning time estimates, the estimated running time for algorithm 1 was calculated using test = tT +((1-P )/P )tp, where tT and tp designate, respectively, the total and parallel times. As in algorithm 2 only one node can "work" at each iteration, the estimated time coincides with the total time. However, we must be aware that it is possible to solve P simultaneous problems using algorithm 2, so the time to solve each problem is, in average, less than the one we obtained (by a factor of P ).

Vertical partition. The simulations of algorithm 3 were done using the network of figure 2, where each node stored one row of a matrix A ∈ R 7×30 (each entry of which was randomly obtained from a Gaussian distribution). To understand how the time estimates were calculated, we must first understand the dynamics of the calculation of new variables. Suppose that, at some iteration t of the nonlinear Gauss-Seidel sweep, node 2 of figure 2 computes x t 2 . The next step consists in sending this variable to the neighbors 1, 3 and 4 which, at iteration t + 1, compute their new variables. After sending them to node 2, this node is again able to compute a new variable, x t+1

2 . Hence, each node computes a new variable every two iterations. This also means that, in average, half of the nodes of a network are "working" (computing new variables) at each iteration. So, the time in parallel can be computed as tp = (2/P )tT , where tT is the total time, whereas the estimated running time is given by the same expression we used for algorithm 1.

General comments. Independently of the complexity of the operations that a node has to perform at each iteration, all the three algorithms require many communications between the nodes before their convergence. Although this depends on several parameters and on the choice of the stopping criteria, we can say that algorithm 1 needs around 700 subgradient iterations, thus 1400P total communications. In the simulations, algorithm 2 required around 3500 total communications between the nodes, and algorithm 3 needed around 8000. In all algorithms, we can increase (resp. reduce) the communication overhead by increasing (resp. decreasing) the number of nodes. Additionally, in algorithm 3 we can explore the con-nectivity of the network: by increasing (resp. decreasing) it, we reduce (resp. increase) significantly the number of iterations of the nonlinear Gauss-Seidel method.

CONCLUSIONS

We propose three algorithms capable of solving the basis pursuit problem (BP) in a distributed environment. That distributivity is characterized by the fact that parts of the matrix A are stored in different nodes, being the nodes able to communicate among themselves.

Two of the algorithms are to be applied in the case where A is partitioned horizontally. The fundamental difference between those two algorithms is in the network topology that they require: one needs a central node and has the advantage that the other (distributed) nodes perform very simple calculations, while the other runs in networks where the nodes are connected circularly. In the latter, the calculations at each node are more complex, but it has the potential to solve P similar problems at the same time.

We also propose an algorithm that is able to solve (BP) when the partition of A is vertical. This algorithm can run in any connected network. This makes this algorithm attractive for sensor networks applications. However, the major disadvantage we found in it was the high number of communications (and thus the time consumed) required to converge, showing that there is margin for improvement.

  Dual problem. Dualizing only the equality constraints in (BBP) leads to the unconstrained dual program maximize L(λ) λ , (DBBP) where L(λ) = b λ + P p=1 Lp(λ) and

Fig. 1 .

 1 Fig. 1. Minimal architecture for the links between the nodes.

Algorithm 1 .

 1 Input: b at the central node (and matrices Ap stored at each node). Output: x at the central node Computation at the central node: 1. Iterate (1) until some stopping criterion is met (a) Send λ k to each node (b) For each node p = 1, . .

3. Send λ to each node 4 .•

 4 For each node p = 1, . . . , P do Computation at node p: Compute the pth fragment of I, say Îp • Return Îp and (Ap) Îp to the central node 5. Receive Îp and (Ap) Îp from each node 6. Concatenate them to obtain Î and A Î 7. Run case 1 or 2 to obtain x

  CBP) where the variable to optimize is x = (x1, . . . , xP ) ∈ (R n ) P and ρ > 0 is a constant. Problem (CBP) is equivalent to (BP) because the graph is connected: thus, x1 = x2 = • • • = xP (we just duplicated the variable). Dual problem. Dualizing only the constraints xi = xj in (CBP) leads to the unconstrained problem maximize L({λ (i,j) }) {λ (i,j) : (i, j) ∈ E} , (DCBP) where L({λ (i,j) }) = inf L({xp}, {λ (i,j) }) : Apxp = bp, for all p and L({xp}, {λ (i,j) }) =

Algorithm 3 .

 3 (Procedure for nodes p = 1, . . . , P ) Input: bp ∈ R mp at each node (and matrices Ap stored at each node). Output: x at each node Iterate on k (method of multipliers) 1. If it receives x t+1 j or x t j from all its neighbors, (a) Form τp = j∈Np (sign(j -p) λ (p,j) -ρxj) (b) Find x t+1 p that solves (12) (using τp) (c) Send x t+1 p to all neighbors (Np)

Table 1 .

 1 Results from the simulations. The numbers shown are an average of 15 random experiments.

		Horizontal partition	Vertical partition
	Feature	Algorithm 1	Algorithm 2	Algorithm 3
	Error in x [%]	2.26	3.69 × 10 -2	6.07 × 10 -1
	Error in f (x) [%]	6.96 × 10 -1 1.96 × 10 -2	7.07 × 10 -2
	Total Time [s]	6.90 × 10 -1 3.66	7.59 × 10 1
	Time in Parallel [s]	5.84 × 10 -1		2.17 × 10 1
	Estimated Time [s]	1.64 × 10 -1 3.66	5.73 × 10 1

It is well known that explicit storage might not be required in some cases, only the ability of implementing the operator Ap and its adjoint A p .

Although this terminology is usually used for the minimization of a convex function, we decided to keep it in our case: the maximization of a concave

function.[START_REF] Boyd | Convex Optimization[END_REF] The vector s is a supergradient (resp. subgradient) of a concave (resp. convex) function f at the point x if f (y) ≤ f (x)+s (y -x) (resp. f (y) ≥ f (x) + s (y -x)) holds for all y.

Partially supported by grants SIPM PTDC/EEA-ACR/73749/2006 and SFRH/BD/33520/2008 (through the Carnegie Mellon/Portugal Program managed by ICTI) from Fundação para a Ciência e Tecnologia; and also by ISR/IST plurianual funding (POSC program, FEDER). This work was also supported by NSF through award 0634967.


Gabriella Bretti
email: gab.bretti@gmail.com

Me Mo

Mat

Massimo Fornasier
email: massimo.fornasier@oeaw.ac.at

Francesca Pitolli
email: pitolli@dmmm.uniroma1.it

Electric current density imaging via an accelerated iterative algorithm with joint sparsity constraints

Many problems in applied sciences require to spatially resolve an unknown electrical current distribution from its external magnetic field. Electric currents emit magnetic fields which can be measured by sophisticated superconducting devices in a noninvasive way. Applications of this technique arise in several fields, such as medical imaging and non-destructive testing, and they involve the solution of an inverse problem.

Assuming that each component of the current density vector possesses the same sparse representation with respect to a preassigned multiscale basis, allows us to apply new regularization techniques to the magnetic inverse problem.

The solution of linear inverse problems with sparsity constraints can be efficiently obtained by iterative algorithms based on gradient steps intertwined with thresholding operations. We test this algorithms to numerically solve the magnetic inverse problem with a joint sparsity constraint.

I. INTRODUCTION

Many problems in applied sciences require to map the spatial distribution of electric currents flowing through a given sample. For example, in neuroscience studies the brain activity can be localized by detecting the regions where neural currents flow. Similarly, in nondestructive testing cracks and corrosion damage in a structure result in perturbations of the current flow.

A given electric current distribution produces a characteristic magnetic field which provides information on the process occurring within the object. Since the magnetic field decreases fast as the distance from the current source to the sensor position increases, the magnetic data can be affected by a high noise. Usually, the measurements are obtained by SQUID magnetometers, that can detect magnetic fields generated by deep currents in a noninvasive way [START_REF] Wikswo | Applications of SQUID magnetometers to biomagnetism and nondestructive evaluation[END_REF].

Of course the magnetic data do not give an immediate image of the current distribution, hence we have to solve the associated inverse problem and regularization techniques are needed [1], [START_REF] Sarvas | Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem[END_REF]. The Tikhonov regularization with quadratic constraint, gives good results when the quantities under observation are equally distributed in time or space [START_REF] Engl | Regularization of inverse problems[END_REF]. However, in the applications we are interested in, the regions where the currents flow are usually small, hence the current distribution is spatially inhomogeneous, and can be represented as a sum of weighted basic currents with only few relevant terms.

To promote a sparse representation regularization techniques based on ℓ 1 -minimization have been introduced, see, e.g., [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF].

In the case the quantity of interest can be modelled as a vector-valued signal with jointly sparse components, then the regularizations introduced in [START_REF] Fornasier | Recovery algorithms for vector valued data with joint sparsity constraints[END_REF], [START_REF] Fornasier | Iterative thresholding algorithms[END_REF] are more suited, see also [START_REF] Bretti | An iterative thresholding algorithm for the neural current imaging[END_REF], [START_REF] Fornasier | Adaptive iterative thresholding algorithms for magnetoencephalography[END_REF], [START_REF] Pitolli | An iterative algorithm with joint sparsity constraints for magnetic tomography[END_REF].

Iterative thresholding algorithms have been proposed for the solution of inverse problems with sparsity constraints, but they have shown rather slow performances in the magnetic tomography problem. In this paper we would like to proposed also for the joint-sparsity constraints a modification inspired by [START_REF] Daubechies | Accelerated projected gradient methods for linear inverse problems with sparsity constraints[END_REF] where the thresholding operations are substituted with suitable projections onto convex sets. These projections turn out to be simply adaptive thresholding steps. For classical sparsity constraints based on ℓ 1 -minimization this approach has shown significant speed improvements in practice.

II. THE FORWARD AND INVERSE MAGNETIC PROBLEMS

From the quasi-static Maxwell's equations it follows that the magnetic field B generated by a current density J obeys the Biot-Savart law:

B( J, r) = µ 0 4π V0 J( r ′ ) × ( r -r ′ ) | r -r ′ | 3 d r ′ , (2.1)
where µ 0 is the magnetic permeability in the vacuum and V 0 is the region where the current flows. Usually, the magnetometers measure just the normal component of the magnetic field, whose expression is given by

B e ( J, r) = µ 0 4π V0 e( r) × ( r ′ -r) | r ′ -r| 3 • J( r ′ ) d r ′ , (2.2)
where r is a point outside V 0 and e( r) is the unit normal vector w.r.t. the surface where the magnetometer is located. We remark that the knowledge of the normal component uniquely determines the magnetic field in the outer space [START_REF] Kress | Reconstruction of a current distribution from its magnetic field[END_REF]. Now, let q l , l = 1, . . . , N , be the magnetometer locations and let M = {m 1 , . . . , m N } be the corresponding measurements. In order to reconstruct the current distribution J we would like to minimize the discrepancy ∆( J) :=

G( J) -M 2 R N
, with G( J) = {B e ( J, q 1 ), ..., B e ( J, q N )}.

Unfortunately, there exist silent currents, so that non unique solutions can be expected. Moreover, the magnetic data can be affected by high noise. A regularization mechanism is required both to identify uniquely the solution by taking advantage of possible prior knowledge (in this case the spatial joint-sparsity of the currents) and to remove the noise.

III. THE MAGNETIC INVERSE PROBLEM WITH SPARSITY

CONSTRAINTS

We are interested in applications where the current density is spatially inhomogeneous so that it can be represented as a sum of weighted basic currents with only few significant terms. This means that we can assume J = (J 1 , J 2 , J 3 ) ∈ L 2 (V 0 ; R 3 ) sparsely represented by a suitable dictionary D := (ψ λ ) λ∈Λ , i.e.

J ℓ ≈ λ∈ΛS j ℓ λ ψ λ , j ℓ λ = J ℓ , ψ λ , ℓ = 1, 2, 3, (3.3)
where Λ S ⊂ Λ, the small set of significant coefficients, is the same for all the components. As a dictionary we can choose a stable multiscale basis, for instance a wavelet basis or frame [START_REF] Christensen | An introduction to frames and Riesz bases[END_REF], [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF].

Thus, the magnetic inverse problem with the sparsity constraint Ψ D consists in minimizing the functional

J Ψ ( j, v) = ∆( j) + Ψ D ( j, v), (3.4)
with respect to both j = (j ℓ λ ) λ∈Λ,ℓ=1,2,3 and an auxiliary weight v, restricted to v λ ≥ 0.

We use as sparsity constraint the joint sparsity measure introduced in [START_REF] Fornasier | Recovery algorithms for vector valued data with joint sparsity constraints[END_REF], i.e.

Ψ D ( j, v) = λ∈Λ v λ j λ p + λ∈Λ ω λ j λ 2 2 + λ∈Λ θ λ (ρ λ -v λ ) 2 , (3.5)
where p ≥ 1 and (θ λ ) λ∈Λ , (ρ λ ) λ∈Λ , (ω λ ) λ∈Λ are positive parameter sequences. Here, • p denotes the usual p-norm for vectors in R 3 . In this way, the minimization of J Ψ promotes that all the entries of the vector j λ = (j ℓ λ ) ℓ=1,2,3 have the same sparsity pattern.

By using the decomposition (3.3), the minimum problem reduces to minimization of the functional

J Ψ ( j, v) = T j -M 2 R N + Ψ D ( j, v), (3.6)
where the entries of T j are given by

(T j) l = 3 ℓ=1 λ∈Λ j ℓ λ µ 0 4π V0 e( q l ) × ( r ′ -q l ) | r ′ -q l | 3 ℓ ψ λ ( r ′ ) d r ′ . (3.7)
An efficient iterative algorithm to numerically solve this minimum problem will be illustrated in the following Section.

IV. AN ACCELERATED ITERATIVE THRESHOLDING

ALGORITHM

The minimizer of the functional J Ψ can be approximated by an accelerated projected gradient method, deduced by combining the iterative thresholding algorithm given in [START_REF] Fornasier | Iterative thresholding algorithms[END_REF],

with the accelerated projected gradient method proposed in [START_REF] Daubechies | Accelerated projected gradient methods for linear inverse problems with sparsity constraints[END_REF].

The iterative thresholding algorithm is as follows (cf.

[2], [START_REF] Pitolli | An iterative algorithm with joint sparsity constraints for magnetic tomography[END_REF]).

Algorithm JS

              
Choose the sequences θ, ρ, ω

Pick an arbitrary j (0) ∈ ℓ 2 (Λ; R 3 )

For 0 ≤ k ≤ K do j (k+1) = S (p) θ,ρ,ω j (k) + T * M -T * T j (k)
The operator S

θ,ρ,ω is a vector-valued thresholding operator whose explicit expression depending on the parameters can be found in [START_REF] Fornasier | Iterative thresholding algorithms[END_REF].

Following the approach in [START_REF] Daubechies | Accelerated projected gradient methods for linear inverse problems with sparsity constraints[END_REF] we define the set The projected gradient iteration has the form: where P C θ,ρ,ω (R) is the orthogonal projection onto C θ,ρ,ω (R).

C θ,ρ,ω (R) := { j ∈ ℓ 2 (Λ; R 3 ) : Ψ D ( j, v) ≤ R, for some 0 ≤ v ≤ ρ}. The set C θ,ρ,ω (R) is convex whenever ω λ θ λ ≥ κ 4 for all λ ∈ Λ,
Algorithm PG                Pick an arbitrary j (0) ∈ ℓ 2 (Λ; R 3 ) For 0 ≤ k ≤ K do Choose the descent parameter β (k) j (k+1) = P C θ,ρ,ω (R) j (k) + β (k) (T * M -T * T j (k) )
This projection can be computed in practice by an adaptive thresholding choice: for any j ∈ ℓ 2 (Λ; R 3 ) there exists a η := η( j, θ, ρ, ω, R) > 0 such that

P C θ,ρ,ω (R) ( j) = S (p)
j ′ ,v ′ j -j ′ 2 R 3 + ηΨ D ( j ′ , v) = arg min j ′ ∈C θ,ρ,ω (R),v ′ j -j ′ 2 R 3 + ηΨ D ( j ′ , v) = arg min j ′ ∈C θ,ρ,ω (R) j -j ′ 2 R 3 + ηR = P C θ,ρ,ω (R) ( j).
The fast convergence of Algorithm PG depends on a suitable choice of the adaptive descent parameter: following the arguments in [START_REF] Daubechies | Accelerated projected gradient methods for linear inverse problems with sparsity constraints[END_REF] it is possible to show that the sequence ( j (k) )

weakly converges to the minimizer of J Ψ which lies on the set C θ,ρ,ω (R) if the sequence β (k) satisfies the conditions

         β := sup{β (k) , k ∈ N} < ∞ , inf{β (k) , k ∈ N} ≥ 1 , β (k) T ( j (k+1) -j (k) ) 2 ≤ r j (k+1) -j (k) 2 , k ≥ k 0 , (4.9)
where r := T * T ℓ2→ℓ2 < 1 and k 0 is a suitable index.

Since at each iteration the evaluation of the action of T * T is requested, we wonder whether such evaluation can be compressed and performed in a fast way.

Let M be the matrix whose entries are the coordinates of T * T in the multiscale basis (ψ λ ) λ∈Λ , i.e.

T * T j λ,ℓ = µ∈Λ 3 m=1 j m µ M (λ,ℓ),(µ,m) , (4.10)
for λ ∈ Λ, ℓ = 1, 2, 3. From (3.7) it follows

M (λ,ℓ),(µ,m) := N l=1 (A ℓ,l ψ λ )(A m,l ψ µ ), (4.11) λ, µ ∈ Λ, ℓ, m = 1, 2, 3,
where

A i,l ψ λ = µ 0 4π V0 e( q l ) × ( r ′ -q l ) | r ′ -q l | 3 ℓ ψ λ ( r ′ ) d r ′ . (4.
12)

It can be shown that M has compressibility properties w.r.t. a compactly supported wavelet basis (see [START_REF] Fornasier | Adaptive iterative thresholding algorithms for magnetoencephalography[END_REF]) so that T * T j (k) can be efficiently evaluated.

V. A BIDIMENSIONAL TEST

In the numerical tests the magnetic data are generated by three horizontal bidimensional current dipoles located in the plane Π 0 = {x, y ∈ R, z = 0}. The magnetic field is sampled by 400 magnetometers located on a regular horizontal grid at height δ = 1. Note that we use adimensional measure units in the tests.

This setting can be used to model both the problem of localizing shallow neural sources and non destructive testing of thin structures. In fact, the current dipoles can be viewed both as sources of brain activity and as discontinuities in a given current distribution. We remark that in the bidimensional case the inverse problem has a unique solution, nevertheless, it can be ill-conditioned for the presence of noise.

As a multiscale basis, we choose the Daubechies orthonormal wavelets with d = 4 vanishing moments and discretize the plane Π 0 with 32 pixels for each dimension. Finally, 2 multiscale levels are used for the current decomposition.

We performed some preliminary tests on Algorithm JS with the thresholding parameter tuned following [START_REF] Fornasier | Recovery algorithms for vector valued data with joint sparsity constraints[END_REF], [START_REF] Fornasier | Iterative thresholding algorithms[END_REF]. In particular, for all λ ∈ Λ we choose ρ λ = 10 -4 in the noiseless case, and ρ λ = 4 • 10 -3 in the noisy case. As for the other parameters, we fix θ λ = 10 -4 and ω λ = 0 for all λ ∈ Λ, thus in the sequel we will suppress the indexes θ, ω from the thresholding operator, that will be simply addressed as S (p) ρ . In Fig. 3 the reconstructed current distribution after 100 iterations by soft-thresholding, i.e. Algorithm JS with p = 1, is shown.

In Fig. 5 the current distribution reconstructed is displayed when high white Gaussian noise with linear signal to noise ratio equal to 0.1 is added to the magnetic field (see Fig. 4).

In spite of the high level of noise and of the poor discretization of the plane, 32 linear pixels, it is still possible to localize the current sources.

As a comparison, in Fig. 6 we show localization results obtained by means of quadratic Tikhonov regularization in the case when the data are distorted by high noise. In real fact, quadratic Tikhonov regularization is not able to give a satisfactory localization of the current sources: when the regularization parameter is low (equal to 10 -3 in the example) the current is not reconstructed at all (Fig. 6, left), while a greater regularization parameter blurres the image (Fig. 6, right, where the regularization parameter is chosen by means of the discrepancy principle and is equal to 356).

VI. CONCLUSIONS

The numerical tests presented in the previous Section show that Algorithm JS behaves better than classical quadratic Tikhonov regularization in case of high noise, and represent promising results for the reconstruction of current sources, in spite of the poor discretization with just 32 linear pixels.

We proposed further Algorithm PG as an acceleration of Algorithm JS generalizing the approach introduced in [START_REF] Daubechies | Accelerated projected gradient methods for linear inverse problems with sparsity constraints[END_REF] to the case of a joint-sparsity constraint.

  where κ = 3 for p = 1, and κ = 1 for p ∈ {2, ∞}, compare [9, Proposition 2.1]. See Figure IV for examples of these sets for different parameters, also in case of violation of the convexity.

Fig. 1 .

 1 Fig. 1. The sets C θ,ρ,ω (R) in R 2 (no vector components are considered and Λ contains two indexes only) for different parameters θ, ρ, ω, R.

θη

  ,ηρ,ηω ( j).

( 4 . 8 )

 48 This relationship is promptly seen by using the definition of S (p) θ,ρ,ω : choose η > 0 such that S (p) θ η ,ηρ,ηω ( j) ∈ ∂C θ,ρ,ω (R), i.e., Ψ D (S (p) θ η ,ηρ,ηω ( j), v) = R for some v and Ψ D (S (p) θ η ,ηρ,ηω ( j), v ′ ) ≥ R for all v ′ . Then S (p) θ η ,ηρ,ηω ( j) = arg min

Fig. 2 .

 2 Fig. 2. The magnetic field produced by three current dipoles located in (-0.5, -0.4, 0), (-0.1, -0.4, 0), (0.3, -0.4, 0). The black points represent the magnetometer sites.

Fig. 3 .Fig. 4 .

 34 Fig. 3. The current intensity reconstructed starting from the magnetic field displayed in Fig. 2 by using soft-thresholding with τ = 10 -4 .

Fig. 5 .Fig. 6 .

 56 Fig. 5. The current intensity reconstructed starting from the magnetic field displayed in Fig. 4 by using soft-thresholding with τ = 4 • 10 -3 .