Nicolas Maillet

Guillaume Collet

Thomas Vannier

Dominique Lavenier

Pierre Peterlongo

COMMET: comparing and combining multiple metagenomic datasets

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

NGS revolution enabled the emergence of the metagenomic field where an environment is sequenced instead of an individual or a species, opening the way to a comprehensive understanding of environmental microbial communities. Large metagenomic projects such as MetaSoil [START_REF] Delmont | Structure, fluctuation and magnitude of a natural grassland soil metagenome[END_REF], MetaHit [START_REF] Qin | A human gut microbial gene catalogue established by metagenomic sequencing[END_REF] or Tara Oceans [START_REF] Karsenti | A holistic approach to marine Eco-systems biology[END_REF] witness this evolution. Analyzing of metagenomic data is a major bottleneck. For instance, assembly tests over "simple" simulated metagenomes showed that N50 is only slightly larger than read sizes [START_REF] Pignatelli | Evaluating the fidelity of De Novo short read metagenomic assembly using simulated data[END_REF]. This situation becomes even worth on complex datasets, such as seawater, where millions of distinct species coexist. In this case, biodiversity can be estimated by using statistical approaches [START_REF] Wang | Metacluster 5.0: A two-round binning approach for metagenomic data for low-abundance species in a noisy sample[END_REF] or by mapping reads on reference banks [START_REF] Huson | MEGAN analysis of metagenomic data[END_REF], [START_REF] Markowitz | IMG/M: The integrated metagenome data management and comparative analysis system[END_REF]. Nevertheless, statistical approaches are limited to a few dozens of species with limited differences in their relative abundance. In addition, the mapping approaches are limited to current knowledge contained in reference banks that suffer from their incompleteness and their inherent errors [START_REF] Schnoes | Annotation error in public databases: Misannotation of molecular function in enzyme superfamilies[END_REF].

A key point of substantial metagenomic projects stands in the number of metagenomes they produce. Then, similarities and differences between metagenomes can be exploited as a source of information, measuring external effects like pollution sources, geographic locations, and patient microbial gut environment [START_REF] Qin | A human gut microbial gene catalogue established by metagenomic sequencing[END_REF], [START_REF] Delmont | Describing microbial communities and performing global comparisons in the omic era[END_REF]. A few methods were proposed to compare metagenomes using external information sources such as taxonomic diversity [START_REF] Jaenicke | Comparative and joint analysis of two metagenomic datasets from a biogas fermenter obtained by 454-pyrosequencing[END_REF] or functional content [START_REF] Sommer | Functional characterization of the antibiotic resistance reservoir in the human microflora[END_REF]. However, these methods are biased because as they are based on partial knowledge.

Methods were proposed to compare metagenomes without using any a priori knowledge. These de novo methods use global features like GC content [START_REF] Foerstner | Environments shape the nucleotide composition of genomes[END_REF], genome size [START_REF] Raes | Prediction of effective genome size in metagenomic samples[END_REF] or sequence signatures [START_REF] Jiang | Comparison of metagenomic samples using sequence signatures[END_REF]. Theses methods face limitations as they are based on rough imprecise criteria and as they only compute a similarity distance: they do not extract similar elements between samples. We believe that it is possible to go further by comparing metagenomic samples at the read sequence level. This provides a higher precision distance and, importantly, it provides reads that are similar between datasets or that are specific to a unique dataset, enabling their latter analysis: assembly with better coverage or comparison with other metagenomic samples. Such comparisons may be performed using Blast [START_REF] Altschul | Basic local alignment search tool[END_REF] or Blat [START_REF] Kent | BLAT-The BLAST-Like Alignment Tool[END_REF] like tools. Unfortunately, these methods do not scale up on large comparative metagenomic studies in which hundreds of millions of reads have to be compared to other hundreds of millions of reads. For instance, one can estimate that comparing a hundred of metagenomes each composed by a hundred of millions of reads of size 100 would require centuries of CPU computation. The crAss approach [START_REF] Dutilh | Reference-independent comparative metagenomics using cross-assembly: crAss[END_REF] constructs a reference metagenome by cross assembling reads of all samples. Then, it maps the initial reads on the so obtained contigs and several measures are derived, based on the repartition of mapped reads. This method provides results of high quality. However, due to its assembly and mapping approach, it does not scale up to large metagenomic datasets. Simpler methods such as Triage-Tools [START_REF] Fimereli | TriageTools : tools for partitioning and prioritizing analysis of high-throughput sequencing data[END_REF] or Compareads [START_REF] Maillet | Compareads: comparing huge metagenomic experiments[END_REF] measure the sequence similarity of a read with a databank by counting the number of kmers (words of length k) shared with the databank. Due to memory consumption, TriageTools cannot use k values larger than 15 and is thus limited to small datasets (a few hundred of thousands reads of length 100). The Compareads tool scales up to large datasets with a small memory footprint and acceptable running time. However, applied on large metagenomic projects, this tools generates an important number of large intermediate result files. In practice, applying Compareads to N datasets generates N 2 resulting new datasets, each of the size of the original ones at worse. Additionally, Compareads leads to highly redundant computation raising up the execution time. These drawbacks are serious bottlenecks limiting the practical usage of Compareads.

In this paper, we introduce COMMET ("COmpare Multiple METagenomes"), a fast software that provides a global similarity overview between all datasets of a metagenomic project. COMMET is based on the Compareads philosophy that consists in determining similarity between two metagenomic datasets by extracting common reads using k-mer approach: two reads are considered similar if they share t non-overlapping k-mers (t and k are parameters). A metagenomic project involving N datasets will thus requires the computation of N 2 intersections which is both time-and storage-consuming. To keep computation time as low as possible, the computation of the N 2 intersections has been strongly improved compared to the Compareads approach through an efficient indexing strategy in which each file is fully indexed only once. In addition, to save storage space, intersections between metagenomic datasets are represented as bit vectors. This compact representation reduces the storage space by two orders of magnitude. Moreover, it provides an easy way to filter and sub-sample reads, or to combine various results by applying logical operations. Finally, COMMET computes a clusterization of metagenomic datasets, which is visualized by dendrogram and heatmaps.

II. METHOD

A. Comparing two sets of reads

The COMMET algorithm to compare two sets of read is based on the Compareads [START_REF] Maillet | Compareads: comparing huge metagenomic experiments[END_REF] methodology. It consists in finding reads from a set A that are similar to at least one read from a set B. The similarity between two reads is based on a minimal number t of non-overlapping identical k-mers. This core operation is directed : it provides reads from A similar to reads from B but it does not provides reads from B similar to reads from A. Note that, as explained below, this operation is based on a heuristic. Thus we denote this operation by A ∩ B.

Computing A ∩ B consists in two steps. Firstly, k-mers from B are indexed in a Bloom filter like data-structure [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF]. Secondly, non-overlapping k-mers of reads from A are searched in the Bloom filter. A read r from A sharing t nonoverlapping k-mers with the Bloom filter is considered similar to at least one read from B. However, the algorithm does not check that these k-mers co-occur on a single read from B, which is a source of false positives. Readers are invited to refer to [START_REF] Maillet | Compareads: comparing huge metagenomic experiments[END_REF] for having more details on precision results.

We recall that the following strategy is applied in order to limit the second source of false positives. First A ∩ B is computed. Then, instead of navely computing B ∩ A, B ∩ (A ∩ B) is computed. This limits the indexed reads of A to those already detected as similar to at least one read from B. Finally, the symmetrical operation is performed:

A ∩ (B ∩ (A ∩ B)).
The previously exposed strategy to fully compare sets A and B within three consecutive ∩ operations has also the advantage to limit the indexation effort. Indeed, only the first A ∩ B operation indexes the full set B. The two other operations only index subsets of A and B.

While comparing read samples A and B, the final results of interest are the reads of A similar to reads of B computed by A ∩ (B ∩ (A ∩ B)) and reads of B similar to reads of A computed by B ∩ (A ∩ B). For sake of simplicity, we denotes these two sets as, respectively, A ⋓ B and B ⋓ A.

In the following sections we present the COMMET novelties: represent read subsets with a limited disk space impact, new read filtering and read subsets manipulation features, compare multiple sets of reads, visualize dataset's similarities as heatmaps and dendrogram.

B. Read subsets representation

In COMMET we propose a simple yet compact data structure to represent a read subset: a vector of bits where each bit represents a read of the original read set. This is what we call the "bit vector representation". As shown below, this representation enables to filter and to subsample read files, to represent ∩ (and thus ⋓ ) results and to easily perform logical operation between read subsets.

Note that with such a representation, a bit vector needs hundreds to thousand times less disk space than a classical uncompressed fastq file. Note also that this way of coding read subsets is not limited to the COMMET framework. It may be applied to any other programs that manipulate read subsets. Thus, the COMMET tool includes a C++ library of reusable components to manipulate read subsets.

In the COMMET framework, the bit vector representation is used as inputs and/or outputs of all tools. In particular they are used in the following operations:

1) Read subsampling and filtering: With huge datasets, it may appear necessary to subsample, for instance limiting each read file to a same number m of a few millions reads. This is immediate by creating a bit vector in which only the first m bits are set to 1, while others are set to 0.

Raw NGS reads also usually need to be filtered on several practical characteristics (read size, read complexity, . . . ). Thus, a bit vector is a direct representation of a filtered result: bit values associated to selected reads are set to 1, the others to 0. A combination of subsampling and filtering allows to select only the m first reads that fulfill the filtration criteria.

2) Representing the similar reads: Results of any ∩ operation is represented by a bit vector. Bit values of reads from the query set detected as similar to at least one read from the reference set are set to 1 and the others are set to 0.

3) Compute logical operations on read subsets: The bit vector representation is ideally suited to perform fundamental logical operations. COMMET provides a module to perform the AN D, OR and N OT operations between distinct subsets of a single initial set of reads.

As presented in the simple case study (Section II-F), these operations, although simple, are powerful while dealing with read subsets. They allow to combine comparison results and so to focus on read subsets intersections or exclusions.

These logical operations perform very efficiently, both in terms of execution time and memory footprint. Moreover, it is worth to notice that they do not generate large result files, as results of these logical operations are also represented as bit vectors. This allows to intensively manipulate read subsets with no technical limitations.

C. Dealing with more than two datasets

We recall that the computation of the A ∩ B core operation involves indexation and search. Once the k-mers of the reads from B are indexed, then the k-mers of the reads from A are sequentially search in the index. If more than a threshold number t of such k-mers are find in the index, then the given read from A is considered as similar to a read from B, which means that the associated value in the bit vector is set to 1.

Consider S = {R 1 , . . . , R N } a set of N ≥ 2 read sets. Applying COMMET on the whole S implies that ∀(i, j) ∈ [1, N ] 2 , i < j, three ordered operations are performed:

1) R i ∩ R j 2) R j ⋓ R i = R j ∩ (R i ∩ R j ) 3) R i ⋓ R j = R i ∩ (R j ∩ (R i ∩ R j ))
Note that for each couple (i, j), the order (i, j) or (j, i) only slightly changes the overall results of the three operations. To avoid redundancies, we limit these operations to i < j.

1) Factorizing the indexation:

In practice, applying COM-MET on S implies to perform the R i ∩ R j operations for all i < j. In particular, R 1 ∩ R N . . . R N -1 ∩ R N have to be computed. For these N -1 computations, the k-mer index of R N is the same. To avoid redundancies, the R N index is computed only once and the N -1 remaining sets are compared to R N using this single index. In general, while

R ref (ref ∈ [2, N ]) is indexed, the index is conserved in RAM memory during the computation of the ref -1 comparisons R query ∩ R ref , with query < ref .
2) Results visualization: Comparisons of N ≥ 2 read sets {R 1 , . . . , R N } provide useful metrics that give an overview of the genomic diversity of the studied samples. Those metrics are summarized in three matrices M 1 , M 2 , and M 3 with values calculated as follows:

• M 1 (i, j) = |R i ⋓ R j | • M 2 (i, j) = 100 × |Ri ⋓ Rj | |Ri| • M 3 (i, j) = 100 × |Ri ⋓ Rj |+|Rj ⋓ Ri| |Ri|+|Rj | . M 1 (i, j) with (i, j) ∈ [1, N ] 2
, is the raw number of reads from R i that are similar to at least one read from R j . As read sets may be of different sizes, M 2 (i, j) is the percentage of reads from R i similar to at least one read from R j . Those two first matrices are asymmetrical. M 3 is a symmetrical matrix. M 3 (i, j) is the percentage of similar reads between the two sets with respect to the total number of reads in R i and R j .

For each matrix, a heatmap is generated . Additionally,

D. The COMMET modules

COMMET integrates four independent modules written in C++, all manipulating, as inputs and outputs, the bit vector representation of read subsets. Additionally, COMMET provides a python script (Commet.py) that takes N ≥ 2 read sets, filters them, compares them and generates explicit representations of comparative results, see Section II-E.

1) Filtering and subsampling reads:

Thanks to the first module, f ilter reads, each read of each dataset (fasta or fastq format, gzipped or not) is filtered out according to user-defined criteria: minimal read length, number of undefined bases, and Shannon complexity [START_REF] Shannon | A mathematical theory of communication[END_REF], used to remove low complexity sequences. The result is a bit vector for each input read file. F ilter reads can also subsample each read set by limiting the number of selected reads to a user defined parameter m. The m first reads that passed the filters are selected.

2) Performing the ∩ core operation: The second module, index and search, performs the ∩ core operation, representing results using the bit vector representation. It inputs a set of read sets (the queries) to be searched in an indexed read set (the bank). A read set may be composed of several read files. Each file could be associated to a bit vector. In this latter case, index and search only considers reads whose associated bit values are set to 1.

3) Manipulating read subsets:

The third module, bvop (bit vector operations) inputs one or two bit vectors. In this second case, the two bit vectors should represent subsets of the same initial set. This module performs the N OT operation on a single bit vector, and the AN D, OR, and AN D N OT operations on two bit vectors.

4) From bit vectors to read files:

Given an original read file and its bit vector, the last module, extract reads, generates an explicit representation of any read subset.

E. Automatization for N ≥ 2 read sets

COMMET includes a python script (Commet.py) which inputs N ≥ 2 read sets. This pipeline i) filters reads, given user-defined parameters, ii) compares all-against-all read sets, and iii) outputs a user-friendly visualization of results. The outputs consist in the three matrices in csv format, their heatmaps and a dendrogram as described in Section II-C2. The dendrogram is realized using the hclust R function, computing a hierarchical complete clustering.

F. Combining read subsets use case

III. RESULTS

A. COMMET efficiently compares multiple metagenomes

We tested COMMET on a set of 28 metagenomes from the IMG/M database [START_REF] Markowitz | IMG/M: The integrated metagenome data management and comparative analysis system[END_REF] (see Table I). These 28 metagenomes were compared with options k = 33, t = 2 and m = 10000. Computations were done using COMMET (Commet.py) and Compareads (v1.3.1) on a 2.9 GHz Intel Core i7 processor with 8GB of RAM and a Solid-State Drive. COMMET calculated the 756 intersections in 35 minutes while Compareads took 81 minutes. In this experiment, COMMET is 2x faster than Compareads thanks to its indexing strategy (each file is fully indexed only once). The obtained dendrograms, shown in Figure 2, are biologically coherent. The different fungus samples are grouped together as well as soil samples, marine planktonic communities and insects. The two human gut microbiome samples are far from other species, as well as the two panda gut microbiome samples.

B. Metasoil study

The MetaSoil study focuses on untreated soils of Park Grass Experiment, Rothamsted Research, Hertfordshire, UK. One of the goals of this study is to assess the influence of depth, seasons and extraction procedure on the sequencing [START_REF] Delmont | Accessing the soil metagenome for studies of microbial diversity[END_REF]. To achieve this, the 13 metagenomes from MetaSoil, two other soil metagenomes and a sea water metagenome, were compared at the functional level using MG-RAST [START_REF] Meyer | The metagenomics RAST server -a public resource for the automatic phylogenetic and functional analysis of metagenomes[END_REF]. This approach identified 835 functional subsystems present in at least one of those metagenomes. On Figure 3.a, samples were clustered using the relative number of reads associated with the 835 functions. This figure shows that the extraction procedure correlates with sample clusters: two metagenomic samples processed with the same extraction procedure share more similarities at the functional level than two samples processed with different extraction procedures [START_REF] Delmont | Structure, fluctuation and magnitude of a natural grassland soil metagenome[END_REF].

This study was reproduced with COMMET on all available metagenomes. The generated bit vectors weigh 68MB while Fig. 3.

Dendrograms from MetaSoil study (top, figure from [START_REF] Delmont | Structure, fluctuation and magnitude of a natural grassland soil metagenome[END_REF]) and COMMET analysis (bottom), comparing the 13 MetaSoil samples, an other soil metagenome and a seawater metagenome (Sargasso Sea).

the explicit representation of the fasta results requires 6.4GB. The storage footprint is thus divided by a factor 100. This ratio is even higher if using fastq format or if dealing with larger read files. The COMMET computation time was 828 minutes (the same set treated by Compareads took 2981 minutes).

Although COMMET uses another metric, the produced dendrogram is highly similar to the MetaSoil one (see Fig 3). On both dendrograms, samples coming from direct extraction are clustered together and external metagenomes are far from the MetaSoil's. Moreover, on the COMMET dendrogram, all samples coming from indirect extraction are clustered together, which is not the case in the MetaSoil study. Even if the two comparing methods are different, they lead to the same conclusion: extraction procedures have a critical impact on sequencing.

IV. CONCLUSION

COMMET gives a global similarity overview of all datasets of a large metagenomic project. It performs all-against-all comparisons of N datasets by factorizing indexation phases. Disk I/Os and storage footprint are highly limited thanks to a new read subset representation which reduces the storage space by at least two orders of magnitude compare to explicit fasta or fastq format. Interestingly, this read subset representation is a powerful way to compute extremely fast boolean operations between read subsets without copying large read files. This enables to focus on reads that fulfill several distinct constraints of interest. The advantages of this representation and of the boolean manipulation are not limited to the COMMET framework. Thus, COMMET includes a C++ library of reusable components to manipulate read subsets.

COMMET produces graphical outputs that sum up allagainst-all comparisons results and open the way for further statistical analysis, thanks to the provided similarity matrices.

A future work will consists in quickly identify significant clusters of read sets by applying rougher comparative metrics (such as the GC content) or a statistical framework based on Principal Component Analysis (PCA). Then, COMMET should be used to go further by precisely compute the shared reads between read sets inside clusters, or between clusters.

COMMET is available under the A-GPL license: http://github.com/pierrepeterlongo/commet.

M 3

 3 is used to construct a dendrogram representation by hierarchical clustering (see Fig 2 for an example of a heatmap and a dendrogram generated by COMMET).

Fig. 1 .

 1 Fig. 1. Logical operations on intersections between A, B and C extract read subsets of interest. The blue subset corresponds to the A AN D N OT (A ⋓ B) operation. The green subset corresponds to the (A ⋓ B) AN D (A ⋓ C) operation. The orange subset corresponds to the (B ⋓ A) OR (B ⋓ C) operation. The red subset corresponds to the (C ⋓ A) AN D N OT (C ⋓ B) operation. By using the bvop module, logical operations can be performed between inputs/outputs of the COMMET pipeline output. For instance, reads from A not similar to any read from set B (blue subset of Fig 1) are obtained by first applying N OT (A ⋓ B) operation. Reads from A similar to at least one read from B and one read from C (green subset

of Fig 1 )

 1 , are identified by computing the AN D operation: (A ⋓ B) AN D (A ⋓ C). In the same spirit, reads from B similar to at least one read from A or one read from C (orange subset of Fig 1), are found by computing the OR operation: (B ⋓ A) OR (B ⋓ C). Operations may be combined to obtain more complex results as, for instance, the red subset of Fig 1, representing reads from C similar to at least one read from A, but not similar to any read from B. This would be done by applying the (C ⋓ A) AN D N OT (C ⋓ B) operation.

Fig. 2 .

 2 Fig. 2. Heatmap (left) and dendrogram (right) representation of the results of the comparison of 28 datasets from the IMG/M database. Results are given with t = 2, m = 10000 and k = 33. The heatmap is constructed from the matrix M 2 and is thus asymmetrical. The dendrogram is constructed from the matrix M 3 by the hierarchical clustering procedure available in R (method "complete").

TABLE I .

 I 28 METAGENOMES FROM THE IMG/M DATABASE

	Identifiers	Description
	SWITGRA	Rhizosphere soil from Panicum virgatum
	SUBGIN	Oral TM7 microbial community of Human
	TERMITE2, TERMITE1	Gut microbiome of divers termites
	SOILM, SOILD, SOILL	Soil microbiome from divers locations
	OMIN, MESO, EUPHO	Divers marine planktonic communities
	BEETLE	Dendroctonus ponderosae
	ACOFUNT, ACOFUNB, CLOFUN,	Fungus garden of divers ants
	ACEFUN, TRAFUN, FUNCOMB
	FUNTER	Fungus-growing termite worker
	WALLABY	Forestomach microbiome of tammar wallaby
	RICE	Endophytic microbiome from rice
	SNAIL	Achatina fulica
	SNOCT	Sirex noctilio microbiome
	XALARV, XAAD	Xyleborus affinis microbiome (larvae, adult)
	HGUT7, HGUT8	Human gut community
	PANDA2, PANDA5	Wild panda gut microbiome

ACKNOWLEDGMENT

Authors warmly thank Claire Lemaitre for her precious advices and her help designing the R functions. This work was supported by the french ANR-2010-COSI-004 MAPPI and by the ANR-12-BS02-0008 Colibread projects. Guillaume Collet's work is funded by the investment expenditure program IDEALG 1192 ANR-10-BTBR-02-04-11.


S Becker

M J Fadili

A quasi-Newton proximal splitting method

A new result in convex analysis on the calculation of proximity operators in certain scaled norms is derived. We describe efficient implementations of the proximity calculation for a useful class of functions; the implementations exploit the piece-wise linear nature of the dual problem. The second part of the paper applies the previous result to acceleration of convex minimization problems, and leads to an elegant quasi-Newton method. The optimization method compares favorably against state-of-the-art alternatives. The algorithm has extensive applications including signal processing, sparse recovery and machine learning and classification.

Introduction

Convex optimization has proved to be extremely useful to all quantitative disciplines of science. A common trend in modern science is the increase in size of datasets, which drives the need for more efficient optimization schemes. For large-scale unconstrained smooth convex problems, two classes of methods have seen the most success: limited memory quasi-Newton methods and non-linear conjugate gradient (CG) methods. Both of these methods generally outperform simpler methods, such as gradient descent.

For problems with non-smooth terms and/or constraints, it is possible to generalize gradient descent with proximal gradient descent (which includes projected gradient descent as a sub-cases), which is just the application of the forward-backward algorithm [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF].

Unlike gradient descent, it is not easy to adapt quasi-Newton and CG methods to problems involving constraints and non-smooth terms. Much work has been written on the topic, and approaches generally follow an active-set methodology. In the limit, as the active-set is correctly identified, the methods behave similar to their unconstrained counterparts. These methods have seen success, but are not as efficient or as elegant as the unconstrained versions. In particular, a sub-problem on the active-set must be solved, and the accuracy of this sub-iteration must be tuned with heuristics in order to obtain competitive results.

Problem statement

Let H = (R N , •, • ) equipped with the usual Euclidean scalar product x, y = N i=1 x i y i and associated norm x =

x, x . For a matrix V ∈ R N ×N in the symmetric positive-definite (SDP) cone S ++ (N ), we define H V = (R N , •, • V ) with the scalar product x, y V = x, V y and norm x V corresponding to the metric induced by V . The dual space of H V , under •, • , is H V -1 . We denote I H the identity operator on H.

A real-valued function f : H → R ∪ {+∞} is (0)-coercive if lim x →+∞ f (x) = +∞. The domain of f is defined by dom f = {x ∈ H : f (x) < +∞} and f is proper if dom f = ∅. We say that a real-valued function f is lower semi-continuous (lsc) if lim inf x→x0 f (x) ≥ f (x 0 ). The class of all proper lsc convex functions from H to R ∪ {+∞} is denoted by Γ 0 (H). The conjugate or Legendre-Fenchel transform of f on H is denoted f * .

Our goal is the generic minimization of functions of the form

min x∈H {F (x) f (x) + h(x)} , (P)
where f, h ∈ Γ 0 (H). We also assume the set of minimizers is nonempty (e.g. F is coercive) and that a standard domain qualification holds. We take f ∈ C 1 (R N ) with L-Lipschitz continuous gradient, and we assume h is separable. Write x ⋆ to denote an element of Argmin F (x).

The class we consider covers non-smooth convex optimization problems, including those with convex constraints. Here are some examples in regression, machine learning and classification.

Example 1 (LASSO).

min x∈H 1 2 Ax -b 2 2 + λ x 1 . (1)
Example 2 (Non-negative least-squares (NNLS)).

min x∈H 1 2 Ax -b 2 2 subject to x 0 . ( 2
)
Example 3 (Sparse Support Vector Machines). One would like to find a linear decision function which minimizes the objective

min x∈H 1 m m i=1 L( x, z i + b, y i ) + λ x 1 (3)
where

for i = 1, • • • , m, (z i , y i ) ∈ R N × {±1}
is the training set, and L is a smooth loss function with Lipschitz-continuous gradient such as the squared hinge loss L(ŷ i , y i ) = max(0, 1ŷi y i ) 2 or the logistic loss L(ŷ i , y i ) = log(1 + e -ŷiyi ).

Contributions

This paper introduces a class of scaled norms for which we can compute a proximity operator; these results themselves are significant, for previous results only cover diagonal scaling (the diagonal scaling result is trivial). Then, motivated by the discrepancy between constrained and unconstrained performance, we define a class of limited-memory quasi-Newton methods to solve (P) and that extends naturally and elegantly from the unconstrained to the constrained case. Most well-known quasi-Newton methods for constrained problems, such as L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF], are only applicable to box constraints l ≤ x ≤ u. The power of our approach is that it applies to a wide-variety of useful non-smooth functionals (see §3.1.4 for a list) and that it does not rely on an active-set strategy. The approach uses the zero-memory SR1 algorithm, and we provide evidence that the non-diagonal term provides significant improvements over diagonal Hessians.

2 Quasi-Newton forward-backward splitting

The algorithm

In the following, define the quadratic approximation

Q B k (x) = f (x k ) + ∇f (x k ), x -x k + 1 2 x -x k 2 B , (4)
where B ∈ S ++ (N ).

The standard (non relaxed) version of the forward-backward splitting algorithm (also known as proximal or projected gradient descent) to solve (P) updates to a new iterate x k+1 according to

x k+1 = argmin x Q B k k (x) + h(x) = prox t k h (x k -t k ∇f (x k )) (5)
with

B k = t -1 k I H , t k ∈]0, 2/L[ (typically t k = 1/L unless a line search is used).
Note that this specializes to the gradient descent when h = 0. Therefore, if f is a strictly convex quadratic function and one takes B k = ∇ 2 f (x k ), then we obtain the Newton method. Let's get back to h = 0. It is now well known that fixed B = LI H is usually a poor choice. Since f is smooth and can be approximated by a quadratic, and inspired by quasi-Newton methods, this suggest picking B k as an approximation of the Hessian. Here we propose a diagonal+rank 1 approximation.

Our diagonal+rank 1 quasi-Newton forward-backward splitting algorithm is listed in Algorithm 1 (with details for the quasi-Newton update in Algorithm 2, see §4 for details). These algorithms are listed as simply as possible to emphasize their important components; the actual software used for numerical tests is open-source and available at http://www.greyc.ensicaen.fr/ ˜jfadili/software.html.

s k ← x k -x k-1
3:

y k ← ∇f (x k ) -∇f (x k-1 ) 4:
Compute H k via Algorithm 2, and define

B k = H -1 k . 5:
Compute the rank-1 proximity operator (see §3)

xk+1 ← prox B k h (x k -H k ∇f (x k )) (6) 6:
p k ← xk+1 -x k and terminate if p k < ǫ 7:
Line-search along the ray x k + tp k to determine x k+1 , or choose t = 1. 8: end for

Relation to prior work

First-order methods The algorithm in ( 5) is variously known as proximal descent or iterated shrinkage/thresholding algorithm (IST or ISTA). It has a grounded convergence theory, and also admits over-relaxation factors α ∈ (0, 1) [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

The spectral projected gradient (SPG) [START_REF] Birgin | Nonmonotone spectral projected gradient methods on convex sets[END_REF] method was designed as an extension of the Barzilai-Borwein spectral step-length method to constrained problems. In [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF], it was extended to non-smooth problems by allowing general proximity operators; The Barzilai-Borwein method [START_REF] Barzilai | Two point step size gradient method[END_REF] uses a specific choice of step-length t k motivated by quasi-Newton methods. Numerical evidence suggests the SPG/SpaRSA method is highly effective, although convergence results are not as strong as for ISTA.

FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] is a multi-step accelerated version of ISTA inspired by the work of Nesterov. The stepsize t is chosen in a similar way to ISTA; in our implementation, we tweak the original approach by using a Barzilai-Borwein step size, a standard line search, and restart [START_REF] O'donoghue | Adaptive restart for accelerated gradient schemes[END_REF], since this led to improved performance. Nesterov acceleration can be viewed as an over-relaxed version of ISTA with a specific, non-constant over-relaxation parameter α k .

The above approaches assume B k is a constant diagonal. The general diagonal case was considered in several papers in the 1980s as a simple quasi-Newton method, but never widely adapted. More recent attempts include a static choice B k ≡ B for a primal-dual method [START_REF] Pock | Diagonal preconditioning for first order primal-dual algorithms in convex optimization[END_REF]. A convergence rate analysis of forward-backward splitting with static and variable B k where one of the operators is maximal strongly monotone is given in [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF].

Active set approaches Active set methods take a simple step, such as gradient projection, to identify active variables, and then uses a more advanced quadratic model to solve for the free variables. A well-known such method is L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF][START_REF] Zhu | Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[END_REF] which handles general box-constrained problems; we test an updated version [START_REF] José | Remark on älgorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound constrained optimization ¨[END_REF]. A recent bound-constrained solver is ASA [START_REF] Hager | A new active set algorithm for box constrained optimization[END_REF] which uses a conjugate gradient (CG) solver on the free variables, and shows good results compared to L-BFGS-B, SPG, GENCAN and TRON. We also compare to several active set approaches specialized for ℓ 1 penalties: "Orthant-wise Learning" (OWL) [START_REF] Andrew | Scalable training of l1-regularized log-linear models[END_REF], "Projected Scaled Sub-gradient + Active Set" (PSSas) [START_REF] Schmidt | Fast optimization methods for l1 regularization: A comparative study and two new approaches[END_REF], "Fixed-point continuation + Active Set" (FPC AS) [START_REF] Wen | A fast algorithm for sparse reconstruction based on shrinkage, subspace optimization and continuation[END_REF], and "CG + IST" (CGIST) [START_REF] Goldstein | High-order methods for basis pursuit[END_REF].

Other approaches By transforming the problem into a standard conic programming problem, the generic problem is amenable to interior-point methods (IPM). IPM requires solving a Newton-step equation, so first-order like "Hessian-free" variants of IPM solve the Newton-step approximately, either by approximately solving the equation or by subsampling the Hessian. The main issues are speed and robust stopping criteria for the approximations.

Yet another approach is to include the non-smooth h term in the quadratic approximation. Yu et al. [START_REF] Yu | A quasi-Newton approach to nonsmooth convex optimization problems in machine learning[END_REF] propose a non-smooth modification of BFGS and L-BFGS, and test on problems where h is typically a hinge-loss or related function.

The projected quasi-Newton (PQN) algorithm [START_REF] Schmidt | Optimizing costly functions with simple constraints: A limited-memory projected quasi-Newton algorithm[END_REF][START_REF] Schmidt | Projected Newton-type methods in machine learning[END_REF] is perhaps the most elegant and logical extension of quasi-Newton methods, but it involves solving a sub-iteration. PQN proposes the SPG [START_REF] Birgin | Nonmonotone spectral projected gradient methods on convex sets[END_REF] algorithm for the subproblems, and finds that this is an efficient tradeoff whenever the cost function (which is not involved in the sub-iteration) is relatively much more expensive to evaluate than projecting onto the constraints. Again, the cost of the sub-problem solver (and a suitable stopping criteria for this inner solve) are issues. As discussed in [START_REF] Lee | Proximal Newton-type methods for minimizing convex objective functions in composite form[END_REF], it is possible to generalize PQN to general non-smooth problems whenever the proximity operator is known (since, as mentioned above, it is possible to extend SPG to this case).

Proximity operators and proximal calculus

For space limitation reasons, we only recall essential definitions. More notions, results from convex analysis as well as proofs can be found in the supplementary material. Definition 4 (Proximity operator [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]). Let h ∈ Γ 0 (H). Then, for every x ∈ H, the function z → 1 2 xz 2 + h(z) achieves its infimum at a unique point denoted by prox h x. The uniquelyvalued operator prox h : H → H thus defined is the proximity operator or proximal mapping of h.

Proximal calculus in H

V Throughout, we denote prox V h = (I H V + V -1 ∂h) -1
, where ∂h is the subdifferential of h, the proximity operator of h w.r.t. the norm endowing H V for some V ∈ S ++ (N ). Note that since V ∈ S ++ (N ), the proximity operator prox V h is well-defined. Lemma 5 (Moreau identity in H V ). Let h ∈ Γ 0 (H), then for any x ∈ H

prox V ρh * (x) + ρV -1 • prox V -1 h/ρ •V (x/ρ) = x, ∀ 0 < ρ < +∞ . (7)
Corollary 6.

prox V h (x) = x -V -1 • prox V -1 h * •V (x) . (8
prox V h (x) = D -1/2 • prox h•D -1/2 (D 1/2 x -v) , (9)
where v = αD -1/2 u and α is the unique root of

p(α) = u, x -D -1/2 • prox h•D -1/2 •D 1/2 (x -αD -1 u) + α , (10)
which is a Lipschitz continuous and strictly increasing function on R with Lipschitz constant

1 + i u 2 i /d i . Remark 8.
• Computing prox V h amounts to solving a scalar optimization problem that involves the computation of prox h•D -1/2 . The latter can be much simpler to compute as D is diagonal (beyond the obvious separable case that we will consider shortly). This is typically the case when h is the indicator of the ℓ 1 -ball or the canonical simple. The corresponding projector can be obtained in expected complexity O(N log N ) by simple sorting the absolute values

• It is of course straightforward to compute prox V h * from prox V h either using Theorem 7, or using this theorem together with Corollary 6 and the Sherman-Morrison inversion lemma.

Diagonal+rank-1: Separable case

The following corollary is key to our novel optimization algorithm.

Corollary 9. Assume that h ∈ Γ 0 (H) is separable, i.e. h(x) = N i=1 h i (x i ), and V = D + uu T , where D is diagonal with (strictly) positive diagonal elements d i , and u ∈ R N . Then,

prox V h (x) = prox hi/di (x i -v i /d i ) i , (11)
where v = αu and α is the unique root of

p(α) = u, x -prox hi/di (x i -αu i /d i ) i + α , (12)
which is a Lipschitz continuous and strictly increasing function on R.

Proof .

Proof: Recall that (10) has a unique solution. When prox hi is piecewise affine with k i segments, it is easy to see that p(α) in ( 12) is also piecewise affine with slopes and intercepts changing at the k transition points di ui (x it j )

(i,j)∈{1,...,N }×{1,...,ki}

. To get α ⋆ , it is sufficient to isolate the unique segment that intersects the abscissa axis. This can be achieved by sorting the values of the transition points which can cost in average complexity O(k log k).

Remark 11.

• The above computational cost can be reduced in many situations by exploiting e.g. symmetry of the h ′ i s, identical functions, etc. This turns out to be the case for many functions of interest, e.g. ℓ 1 -norm, indicator of the ℓ ∞ -ball or the positive orthant, and many others; see examples hereafter.

• Corollary 9 can be extended to the "block" separable (i.e. separable in subsets of coordinates) when D is piecewise constant along the same block indices.

Semi-smooth Newton method

In many situations (see examples below), the root of p(α) can be found exactly in polynomial complexity. If no closed-form is available, one can appeal to some efficient iterative method to solve (10) (or ( 12)). As p is Lipschitz-continuous, hence so-called Newton (slantly) differentiable, semi-smooth Newton are good such solvers, with the proviso that one can design a simple slanting function which can be algorithmically exploited.

The semi-smooth Newton method for the solution of ( 10) can be stated as the iteration

α t+1 = α t -g(α t ) -1 p(α t ) , (13)
where g is a generalized derivative of p. Proposition 12 (Generalized derivative of p). If prox h•D -1/2 is Newton differentiable with generalized derivative G, then so is the mapping p with a generalized derivative

g(α) = 1 + u, D -1/2 • G(D 1/2 x -αD -1/2 u) • D -1/2 u .
Furthermore, g is nonsingular with a uniformly bounded inverse on R. Thus, as p is Newton differentiable with nonsingular generalized derivative whose inverse is also bounded, the general semi-smooth Newton convergence theorem implies that (13) converges superlinearly to the unique root of (10).

Examples

Many functions can be handled very efficiently using our results above. For instance, Table 1 summarizes a few of them where we can obtain either an exact answer by sorting when possible, or else by minimizing w.r.t. to a scalar variable (i.e. finding the unique root of ( 10)).

A primal rank 1 SR1 algorithm

Following the conventional quasi-Newton notation, we let B denote an approximation to the Hessian of f and H denote an approximation to the inverse Hessian. All quasi-Newton methods update an approximation to the (inverse) Hessian that satisfies the secant condition:

H k y k = s k , y k = ∇f (x k ) -∇f (x k-1 ), s k = x k -x k-1 (14)
Algorithm 1 follows the SR1 method [START_REF] Broyden | Quasi-Newton methods and their application to function minimization[END_REF], which uses a rank-1 update to the inverse Hessian approximation at every step. The SR1 method is perhaps less well-known than BFGS, but it has the crucial property that updates are rank-1, rather than rank-2, and it is described "[SR1] has now taken its place alongside the BFGS method as the pre-eminent updating formula." [START_REF] Gould | Seminal papers in nonlinear optimization[END_REF].

We propose two important modifications to SR1. The first is to use limited-memory, as is commonly done with BFGS. In particular, we use zero-memory, which means that at every iteration, a new diagonal plus rank-one matrix is formed. The other modification is to extend the SR1 method to the general setting of minimizing f + h where f is smooth but h need not be smooth; this further generalizes the case when h is an indicator function of a convex set. Every step of the algorithm replaces f with a quadratic approximation, and keeps h unchanged. Because h is left unchanged, the subgradient of h is used in an implicit manner, in comparison to methods such as [START_REF] Yu | A quasi-Newton approach to nonsmooth convex optimization problems in machine learning[END_REF] that use an approximation to h as well and therefore take an explicit subgradient step.

Choosing H 0 In our experience, the choice of H 0 is best if scaled with a Barzilai-Borwein spectral step length τ BB2 = s k , y k / y k , y k (15) (we call it τ BB2 to distinguish it from the other Barzilai-Borwein step size τ BB1 = s k , s k / s k , y k τ BB2 ).

In SR1 methods, the quantity s k -H 0 y k , y k must be positive in order to have a well-defined update for u k . The update is:

H k = H 0 + u k u T k , u k = (s k -H 0 y k )/ s k -H 0 y k , y k . ( 16
)
Algorithm 2: Sub-routine to compute the approximate inverse Hessian H k Require: k, s k , y k , 0 < γ < 1, 0 < τ min < τ max 1: if k = 1 then 2:

H 0 ← τ I H where τ > 0 is arbitrary 3:

u k ← 0 4: else 5: τ BB2 ← s k ,y k
u k ← (s k -H 0 y k )/ s k -H 0 y k , y k ).
12:

end if 13: end if 14: return H k = H 0 + u k u T k {B k = H -1
k can be computed via the Sherman-Morrison formula}

For this reason, we choose H 0 = γτ BB2 I H with 0 < γ < 1, and thus 0 ≤ s k -H 0 y k , y k = (1γ) s k , y k . If s k , y k = 0, then there is no symmetric rank-one update that satisfies the secant condition. The inequality s k , y k > 0 is the curvature condition, and it is guaranteed for all strictly convex objectives. Following the recommendation in [START_REF] Nocedal | Numerical Optimization[END_REF], we skip updates whenever s k , y k cannot be guaranteed to be non-zero given standard floating-point precision.

A value of γ = 0.8 works well in most situations. We have tested picking γ adaptively, as well as trying H 0 to be non-constant on the diagonal, but found no consistent improvements.

Numerical experiments and comparisons

1 2 Ax + -Ax --b 2 + λ1 T (x + + x -)
and then x is recovered via x = x +x -. With such a formulation solvers such as L-BFGS-B are applicable. However, this constrained problem has twice the number of variables, and the Hessian of the quadratic part changes from A T A to à = A T A -A T A -A T A A T A which necessarily has (at least) n degenerate 0 eigenvalues and adversely affects solvers.

A similar situation occurs with the hinge-loss function. Consider the shifted and reversed hinge loss function h(x) = max(0, x). Then one can split x = x +x -, add constraints x + 0, x - 0, and replace h(x) with 1 T (x + ). As before, the Hessian gains n degenerate eigenvalues.

We compared our proposed algorithm on the LASSO problem. The first example, in Fig. 1a, is a typical example from compressed sensing that takes A ∈ R m×n to have iid N (0, 1) entries with m = 1500 and n = 3000. We set λ = 0.1. L-BFGS-B does very well, followed closely by our proposed SR1 algorithm and PSSas. Note that L-BFGS-B and ASA are in Fortran and C, respectively (the other algorithms are in Matlab).

Our second example uses a square operator A with dimensions n = 13 3 = 2197 chosen as a 3D discrete differential operator. This example stems from a numerical analysis problem to solve a discretized PDE as suggested by [START_REF] Fletcher | On the Barzilai-Borwein method[END_REF]. For this example, we set λ = 1. For all the solvers, we use the same parameters as in the previous example. Unlike the previous example, Fig. 1b now shows that L-BFGS-B is very slow on this problem. The FPC-AS method, very slow on the earlier test, is now the fastest. However, just as before, our SR1 method is nearly as good as the best algorithm. This robustness is one benefit of our approach, since the method does not rely on active-set identifying parameters and inner iteration tolerances.

Conclusions

In this paper, we proposed a novel variable metric (quasi-Newton) forward-backward splitting algorithm, designed to efficiently solve non-smooth convex problems structured as the sum of a smooth term and a non-smooth one. We introduced a class of weighted norms induced by a diagonal+rank 1 symmetric positive definite matrices, and proposed a whole framework to compute a proximity operator in the weighted norm. The latter result is distinctly new and is of independent interest. We also provided clear evidence that the non-diagonal term provides significant acceleration over diagonal matrices.

The proposed method can be extended in several ways. Although we focused on forward-backward splitting, our approach can be easily extended to the new generalized forward-backward algorithm of [START_REF] Raguet | Generalized forward-backward splitting[END_REF]. However, if we switch to a primal-dual setting, which is desirable because it can handle more complicated objective functionals, updating B k is non-obvious. Though one can think of non-diagonal pre-conditioning methods.

Another improvement would be to derive efficient calculation for rank-2 proximity terms, thus allowing a 0-memory BFGS method. We are able to extend (result not presented here) Theorem 7 to diagonal+rank r matrices. However, in general, one must solve an r-dimensional inner problem using the semismooth Newton method.

A final possible extension is to take B k to be diagonal plus rank-1 on diagonal blocks, since if h is separable, this is still can be solved by our algorithm (see Remark 10). The challenge here is adapting this to a robust quasi-Newton update. For some matrices that are well-approximated by low-rank blocks, such as H-matrices [START_REF] Hackbusch | A sparse matrix arithmetic based on H-matrices. Part I: Introduction to H-matrices[END_REF], it may be possible to choose B k ≡ B to be a fixed preconditioner.

) 3 . 1 . 1

 311 Diagonal+rank-1: General case Theorem 7 (Proximity operator in H V ). Let h ∈ Γ 0 (H) and V = D + uu T , where D is diagonal with (strictly) positive diagonal elements d i , and u ∈ R N . Then,

y k 2 {Barzilai-Borwein step length} 6 :

 26 Project τ BB2 onto [τ min , τ max ]

7 : 8 :

 78 H 0 ← γτ BB2 I H if s k -H 0 y k , y k ≤ 10 -8 y k 2 s k -H 0 y k 2 then

Figure 1 :

 1 Figure 1: (a) is first LASSO test, (b) is second LASSO test Consider the unconstrained LASSO problem (1). Many codes, such as [27] and L-BFGS-B [2], handle only non-negativity or box-constraints. Using the standard change of variables by introducing the positive and negative parts of x, the LASSO can be recast as min x+,x-0

Table 1 :

 1 Summary of functions which have efficiently computable rank-1 proximity operatorsProof:This follows from linearity and the chain rule[START_REF] Griesse | A semismooth Newton method for Tikhonov functionals with sparsity constraints[END_REF] Lemma 3.5]. The second statement follows strict increasing monotonicity of p as established in Theorem 7.

	Function h	Algorithm
	ℓ 1 -norm	Separable: exact in O(N log N )
	Hinge	Separable: exact in O(N log N )
	ℓ ∞ -ball	Separable: exact in O(N log N ) from ℓ 1 -norm by Moreau-identity
	Box constraint	Separable: exact in O(N log N )
	Positivity constraint Separable: exact in O(N log N )
	ℓ 1 -ball ℓ ∞ -norm	Nonseparable: semismooth Newton and prox h•D -1/2 costs O(N log N ) Nonseparable: from projector on the ℓ 1 -ball by Moreau-identity
	Canonical simplex max function	Nonseparable: semismooth Newton and prox h•D -1/2 costs O(N log N ) Nonseparable: from projector on the simplex by Moreau-identity

Acknowledgments

SB would like to acknowledge the Fondation Sciences Mathématiques de Paris for his fellowship.


François-Xavier Dupé
email: francois-xavier.dupe@cea.fr

Jalal M Fadili
email: jalal.fadili@greyc.ensicaen.fr

Jean-Luc Starck
email: jean-luc.starck@cea.fr

Linear inverse problems with various noise models and mixed regularizations

Keywords: Inverse Problems, Poisson noise, Gaussian noise, Multiplicative noise, Duality, Proximity operator, Sparsity

In this paper, we propose two algorithms to solve a large class of linear inverse problems when the observations are corrupted by various types of noises. A proper data fidelity term (log-likelihood) is introduced to reflect the statistics of the noise (e.g. Gaussian, Poisson, multiplicative, etc.) independently of the degradation operator. On the other hand, the regularization is constructed through different terms reflecting a priori knowledge on the images. Piecing together the data fidelity and the prior terms, the solution to the inverse problem is cast as the minimization of a composite non-smooth convex objective functional. We establish the well-posedness of the optimization problem, characterize the corresponding minimizers for different kind of noises. Then we solve it by means of primal and primal-dual proximal splitting algorithms originating from the field of non-smooth convex optimization theory. Experimental results on deconvolution, inpainting and denoising with some comparison to prior methods are also reported.

INTRODUCTION

Previous work A lot of works have already been dedicated to linear inverse problems with Gaussian noise (see [START_REF] Starck | Sparse Signal and Image Processing: Wavelets, Curvelets and Morphological Diversity[END_REF] for a comprehensive review), while linear inverse problems in presence of other kind of noise such as Poisson noise have attracted less interest, presumably because noises properties are more complicated to handle. Such inverse problems have however important applications in imaging such as restoration (e.g. deconvolution in medical and astronomical imaging), or reconstruction (e.g. computerized tomography).

Since the work for Gaussian noise by [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraints[END_REF], many other methods have appeared for managing linear inverse problem with sparsity regularization. But they limited to the Gaussian case. In the context of Poisson linear inverse problems using sparsity-promoting regularization, a few recent algorithms have been proposed. For example, [START_REF] Dupé | A proximal iteration for deconvolving poisson noisy images using sparse representations[END_REF] stabilize the noise and proposed a family of nested schemes relying upon proximal splitting algorithms (Forward-Backward and Douglas-Rachford) to solve the corresponding optimization problem. The work of [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] is in the same vein. These methods may be extended to other kind of noise. However, nested algorithms are time-consuming since they necessitate to sub-iterate. Using the augmented Lagrangian method with the alternating method of multipliers algorithm (ADMM), which is nothing but the Douglas-Rachford splitting applied to the Fenchel-Rockafellar dual problem, [START_REF] Figueiredo | Restoration of poissonian images using alternating direction optimization[END_REF] presented a deconvolution algorithm with TV and sparsity regularization, and [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF] a denoising algorithm for multiplicative noise. This scheme however necessitates to solve a least-square problem which can be done explicitly only in some cases.

Contributions In this paper, we propose a framework for solving linear inverse problems when the observations are corrupted by various types of noise. In order to form the data fidelity term, we take the exact likelihood associated to the noise model. As a prior, the images are assumed to comply with several regularity properties or/and constraints reflecting knowledge about the original image. The solution to the inverse problem is cast as the minimization of a composite non-smooth convex functional, for which we prove well-posedness of the optimization problem, characterize the corresponding minimizers, and solve them by means of primal and primal-dual proximal splitting algorithms originating from the realm of non-smooth convex optimization theory. Convergence of the algorithms is also shown. Experimental results and comparison to other algorithms on deconvolution are finally conducted.

Notation and terminology

Let H a real Hilbert space, here a finite dimensional vector subspace of R n . We denote by . the norm associated with the inner product in H, and I is the identity operator on H. . p , p ≥ 1 is the ℓp norm. x and α are respectively reordered vectors of image samples and transform coefficients. We denote by ri C the relative interior of a convex set Let x ∈ H be a n-pixels image. x can be written as the superposition of elementary atoms ϕγ parameterized by γ ∈ I such that x = P γ∈I αγϕγ = Φα, |I| = L, L n. We denote by Φ : H ′ → H the dictionary (typically a frame of H), whose columns are the atoms all normalized to a unit ℓ2-norm

C. A real-valued function f is coercive, if lim x →+∞ f (x) = +∞, and is proper if its domain is non-empty dom f = {x ∈ H | f (x) < +∞} = ∅. Γ0(H)

PROBLEM STATEMENT

Consider the image formation model where an input image of n pixels x is indirectly observed through the action of a bounded linear operator H : H → K, with K a real Hilbert space (usually a subspace of R m , m > 0), and contaminated by a noise ε through a composition operator ⊙ (e.g. addition),

y ∼ Hx ⊙ ε . (1)
The linear inverse problem at hand is to reconstruct x from the observed image y.

A natural way to attack this problem would be to adopt a maximum a posteriori (MAP) bayesian framework with an appropriate likelihood function -the distribution of the observed data y given an original x -reflecting the statistics of the noise. As a prior, the image is supposed to be economically (sparsely) represented in a pre-chosen dictionary Φ as measured by a sparsity-promoting penalty Ψ supposed throughout to be convex but non-smooth, e.g. the ℓ1 norm.

Gaussian noise case

For Gaussian noise, we consider the following formation model,

y = Hx + ε , (2)
where ε ∼ N (0, σ 2 ).

>From the probability density function, the negative loglikelihood writes:

fGaussian : η ∈ H → η -y 2 2 /(2σ 2 ) . (3)
>From this function, we can directly derive the following result, Proposition 1. fGaussian is a proper, strictly convex and lsc function.

Poisson noise case

The observed image is then a discrete collection of counts y = (y[i]) 1 i n which are bounded, i.e. y ∈ ℓ∞. Each count y[i] is a realization of an independent Poisson random variable with a mean (Hx)i. Formally, this writes in a vector form as

y ∼ P(Hx) . ( 4
)
>From the probability density function of a Poisson random variable, the likelihood writes:

p(y|x) = Y i ((Hx)[i]) y[i] exp (-(Hx)[i]) y[i]! . (5)
Taking the negative log-likelihood, we arrive at the following data fidelity term:

fPoisson : η ∈ H → n X i=1 fp(η[i]), (6)
if

y[i] > 0, fp(η[i]) =  -y[i] log(η[i]) + η[i] if η[i] > 0, +∞ otherwise, if y[i] = 0, fp(η[i]) =  η[i] if η[i] ∈ [0, +∞), +∞ otherwise.
Using classical results from convex theory, we can show that, Proposition 2. fPoisson is a proper, convex and lsc function. fPoisson is strictly convex if and only if ∀i ∈ {1, . . . , n}, y[i] = 0.

Multiplicative noise

We consider here the same context of multiplicative as in [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF]. With multiplicative noise, a usual approach to improve the signal to noise ratio consists in averaging independent observation of the same resolution. When considering SAR/SAS system, this method is called multilook, i.e. M -look in the case of the averaging of M images. For fully developed speckle, the averaged images are Gamma distributed,

y = xε, ε ∼ Γ(M, 1/M ) . (7)
In order to simplify the problem, the logarithm of the observation is considered,

log(y) = log(x) + log(ε) = z + ω . (8)
And in [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF], the authors proof that the anti log-likelihood yields,

f Multi : η ∈ H → M n X i=1 (z[i] + exp(log(y[i]) -z[i]) . (9)
Using classical results from convex theory, we can directly derive, Proposition 3. f Multi is a proper, strictly convex and lsc function.

Optimization problem

Our aim is then to solve the following optimization problems, under a synthesis-type sparsity prior 1 where H ′ is the Hilbert space defined by the dictionary, argmin

α∈H ′ J(α), J : α → f1 • H • Φ(α) + K X i=1 Ri(α) . (P f 1 ,γ,ψ )
The data fidelity term f1 reflect the noise statistics and the Ri, 1 i K the K prior terms. In this paper, we restrict to the case K = 2, with K1 = Ψ, the penalty function Ψ :

α → P L i=0 ψi(α[i]
) which is positive, additive, and chosen 1 Our framework and algorithms extend to an analysis-type prior just as well.

to enforce sparsity, γ > 0 is a regularization parameter and K2 = ıC • Φ the indicator function of the convex set C (e.g. the positive orthant for Poissonian data).

For the rest of the paper, we assume that f1 is a proper, convex and lsc function, i.e. f1 ∈ Γ0(H). This is true for many kind of noises including Poisson, Gaussian, Laplacian. . . (see [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] for others examples). >From the objective in (P f 1 ,γ,ψ ), we get the following, Proposition 4. (i) f1 is a convex functions and so are f1•H and f1

•H•Φ. (ii) Suppose that f1 is strictly convex, then f1 • H • Φ re- mains strictly convex if Φ is an orthobasis and ker(H) = ∅. (iii) Suppose that (0, +∞) ∩ H ([0, +∞)) = ∅. Then J ∈ Γ0(H).
2.5 Well-posedness of (P f 1 ,γ,ψ )

Let M be the set of minimizers of problem (P f 1 ,γ,ψ ). Suppose that Ψ is coercive. Thus J is coercive. Therefore, the following holds:

Proposition 5.

(i) Existence: (P f 1 ,γ,ψ ) has at least one solution, i.e. M = ∅. (ii) Uniqueness: (P f 1 ,γ,ψ ) has a unique solution if Ψ is strictly convex, or under (ii) of Proposition 4.

ITERATIVE MINIMIZATION ALGORITHMS 3.1 Proximal calculus

We are now ready to describe the proximal splitting algorithms to solve (P f 1 ,γ,ψ ). At the heart of the splitting framework is the notion of proximity operator. Definition 6 ( [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]). Let F ∈ Γ0(H). Then, for every x ∈ H, the function y → F (y) + xy 2 2 /2 achieves its infimum at a unique point denoted by prox F x. The operator prox F : H → H thus defined is the proximity operator of F .

Then, the proximity operator of the indicator function of a convex set is merely its orthogonal projector. One important property of this operator is the separability property:

Lemma 7 ([7]). Let F k ∈ Γ0(H), k ∈ {1, • • • , K} and let G : (x k ) 1 k K → P k F k (x k ). Then prox G = (prox F k ) 1 k K .
For Gaussian noise, we can easily prove that with f1 as defined in (3), Lemma 8. Let y be the observation, the proximity operator associated to fGaussian (i.e. the Gaussian anti loglikelihood) is,

prox βf Gaussian x = βy + σ 2 x β + σ 2 . ( 10
)
The following result can be proved easily by solving the proximal optimization problem in Definition 6 with f1 as defined in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], see also [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF].

Lemma 9. Let y be the count map (i.e. the observations), the proximity operator associated to fPoisson (i.e. the Poisson anti log-likelihood) is,

prox βf Poisson x = x[i] -β + p (x[i] -β) 2 + 4βy[i] 2 ! 1 i n . (11)
As with multiplicative noise f Multi involves the exponential, we need the W-Lambert function [START_REF] Corless | On the Lambert W function[END_REF] in order to derive a closed form of the proximity operator, Lemma 10. Let y be the observations, the proximity operator associated to f Multi is,

prox βf Multi x = `x[i] -βM - W (-βM exp(x[i] -log(y[i]) -βM )) ´1 i n , ( 12
)
where W is the W-Lambert function, i.e. the function such that W(a) exp(W(a)) = a, a ∈ R.

We now turn to prox γΨ which is given by Lemma 7 and the following result:

Theorem 11 ( [START_REF] Fadili | Inpainting and zooming using sparse representations[END_REF]). Suppose that ∀ i: (i) ψi is convex even-symmetric, non-negative and non-decreasing on R + , and ψi(0) = 0; (ii) ψi is twice differentiable on R \ {0}; (iii) ψi is continuous on R, and admits a positive right derivative at zero ψ

′ i + (0) = lim h→0 + ψ i (h) h > 0.
Then, the proximity operator prox δψ i (β) = α(β) has exactly one continuous solution decoupled in each coordinate β[i] :

α[i] = ( 0 if |β[i]| δψ ′ i + (0) βi -δψ ′ i ( α[i]) if |β[i]| > δψ ′ i + (0) (13)
Among the most popular penalty functions ψi satisfying the above requirements, we have ψi(α[i]) = |α[i]| , ∀ i, in which case the associated proximity operator is soft-thresholding, denoted ST in the sequel.

Splitting on the primal problem

Splitting for sums of convex functions

Suppose that the objective to be minimized can be expressed as the sum of K functions in Γ0(H), verifying domain qualification conditions:

argmin x∈H F (x) = K X k=1 F k (x) ! . ( 14
)
Proximal splitting methods for solving [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF] are iterative algorithms which may evaluate the individual proximity operators prox F k , supposed to have an explicit convenient structure, but never proximity operators of sums of the F k .

Splitting algorithms have an extensive literature since the 1970's, where the case K = 2 predominates. Usually, splitting algorithms handling K > 2 have either explicitly or implicitly relied on reduction of (18) to the case K = 2 in the product space H K . For instance, applying the Douglas-Rachford splitting to the reduced form produces Spingarn's method, which performs independent proximal steps on each F k , and then computes the next iterate by essentially averaging the individual proximity operators. The scheme described in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] is very similar in spirit to Spingarn's method, with some refinements.

Application to noisy inverse problems

Problem (P f 1 ,γ,ψ ) is amenable to the form [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF], by wisely introducing auxiliary variables. As (P f 1 ,γ,ψ ) involves two linear operators (Φ and H), we need two of them, that we define as x1 = Φα and x2 = Hx1. The idea is to get rid of the composition of Φ and H. Let the two linear operators L1 = [I 0 -Φ] and L2 = [-H I 0]. Then, the optimization problem (P f 1 ,γ,ψ ) can be equivalently written:

argmin (x 1 ,x 2 ,α)∈H×K×H ′ f1(x2) + ıC(x1) + γΨ(α) | {z } G(x 1 ,x 2 ,α) + ( 15
)
ı ker L 1 (x1, x2, α) + ı ker L 2 (x1, x2, α) , (16)
where ı ker L i are the indicatrice function of the kernel space of the operator Li, in others words these two indicatrice functions will enforce the equality constraints x1 = Φα and x2 = Hx1. Notice that in our case K = 3 by virtue of separability of the proximity operator of G in x1, x2 and α; see Lemma 7.

Algorithm 1: Primal scheme for solving (P f 1 ,γ,ψ ). Parameters: The observed image y, the dictionary Φ, number of iterations Niter, µ > 0 and regularization parameter γ > 0.

Initialization: ∀i ∈ {1, 2, 3}, p (0,i) = (0, 0, 0) T . z0 = (0, 0, 0) T .

Main iteration:

For t = 0 to Niter -1,

• Data fidelity (Lemmas 8, 9 and 10): [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF]). • Sparsity-penalty (Lemma 11): ξ (t,1) [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] = prox µγΨ/3 (p (t,1) [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]). • Positivity constraint: ξ (t,1) [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] = PC(p (t,1) [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]).

ξ (t,1) [1] = prox µf 1 /3 (p (t,1)
• Auxiliary constraints with L1 and L2: (Lemma 12):

ξ (t,2) = P ker L 1 (p (t,2) ), ξ (t,3) = P ker L 2 (p (t,3) ). • Average the proximity operators: ξt = (ξ (t,1) + ξ (t,2) + ξ (t,3) )/3. • Choose θt ∈]0, 2[. • Update the components: ∀i ∈ {1, 2, 3}, p (t+1,i) = p (t,i) + θt(2ξt -zt -ξ (t,i) ). • Update the coefficients estimate: zt+1 = zt + θt(ξt -zt).
End main iteration Output:

Reconstructed image x ⋆ = zN iter [0].
The proximity operators of f1 and Ψ are easily accessible through Lemmas 8, 9, 10 and 11. The projector onto C is trivial for most of the case (e.g. positive orthant, closed interval). It remains now to compute the projector on ker Li, i = 1, 2, which by well-known linear algebra arguments, is obtained from the projector onto the image of L * i .

Lemma 12. The proximity operator associated to ı ker L i is

P ker L i = I -L * i (Li • L * i ) -1 Li . ( 17
)
The inverse in the expression of P ker L 1 is (I+Φ•Φ T ) -1 can be computed efficiently when Φ is a tight frame. Similarly, for L2, the inverse writes (I+H•H * ) -1 , and its computation can be done in the domain where H is diagonal; e.g. Fourier for convolution or pixel domain for mask.

Finally, the main steps of our primal scheme are summarized in Algorithm 1. Its convergence is a corollary of [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF][Theorem 3.4].

Proposition 13. Let (zt) t∈N be a sequence generated by Algorithm 1. Suppose that Proposition 4-(iii) is verified, and P t∈N θt(2θt) = +∞. Then (zt) t∈N converges to a (nonstrict) global minimizer of (P f 1 ,γ,ψ ).

Splitting on the dual: Primal-dual algorithm

Our problem (P f 1 ,γ,ψ ) can also be rewritten in the form,

argmin α∈H ′ F • K(α) + γΨ(α) (18)
where now

K = " H • Φ Φ « and F : (x1, x2) → f1(x1) + ıC(x2
). Again, one may notice that the proximity operator of F can be directly computed using the separability in x1 and x2.

Recently, a primal-dual scheme, which turns to be a preconditioned version of ADMM, to minimize objectives of the form (18) was proposed in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. Transposed to our setting, this scheme gives the steps summarized in Algorithm 2.

Adapting the arguments of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], convergence of the sequence (αt) t∈N generated by Algorithm 2 is ensured.

Proposition 14. Suppose that Proposition 4-(iii) holds. Let ζ = |||Φ||| 2 (1 + |||H||| 2 ), choose τ > 0 and σ such that στ ζ < 1, and let (αt) t∈R as defined by Algorithm 2. Then, (α) t∈N converges to a (non-strict) global minimizer (P f 1 ,γ,ψ ) at the rate O(1/t) on the restricted duality gap.

Discussion

Algorithm 1 and 2 share some similarities, but exhibit also important differences. For instance, the primal-dual algorithm enjoys a convergence rate that is not known for the primal algorithm. Furthermore, the latter necessitates two operator inversions that can only be done efficiently for some Φ and H, while the former involves only application of these linear operators and their adjoints. Consequently, Algorithm 2 can virtually handle any inverse problem with a bounded linear H. In case where the inverses can be done efficiently, e.g. deconvolution with a tight frame, both algorithms have comparable computational burden. In general, if other regularizations/constraints are imposed on the solution, in the form of additional proper lsc convex terms that would appear in (P f 1 ,γ,ψ ), both algorithms still apply by introducing wisely chosen auxiliary variables. ´.

• Update the coefficients estimate: ᾱt+1 = 2αt+1αt

End main iteration Output: Reconstructed image x ⋆ = ΦαN iter .

EXPERIMENTAL RESULTS

Deconvolution under Poisson noise

Our algorithms were applied to deconvolution. In all experiments, Ψ was the ℓ1-norm. Table 1 summarizes the mean absolute error (MAE) and the execution times for an astronomical image, where the dictionary consisted of the wavelet transform (here the translation-invariant discret wavelet transform which form a tight frame and so is a redundant transform) and the PSF was that of the Hubble telescope. Our algorithms were compared to state-of-theart alternatives in the literature. In summary, flexibility of our framework and the fact that Poisson noise was handled properly, demonstrate the capabilities of our approach, and allow our algorithms to compare very favorably with other competitors. The computational burden of our approaches is also among the lowest, typically faster than the PIDAL algorithm [START_REF] Figueiredo | Restoration of poissonian images using alternating direction optimization[END_REF]. Fig. 1 displays the objective as a function of the iteration number and times (in s). We can clearly see that Algorithm 2 converges faster than Algorithm 1.

RL-MRS [

Inpainting with Gaussian noise

We also applied our algorithms to inpainting with Gaussian noise. In all experiments, Ψ was the ℓ1-norm. Fig 2 summarizes the results with the PSNR and the execution times for the Cameraman, where the dictionary consisted of the same wavelet transform as used for deconvolution and the mask was create from a random process (here with about 34% of missing pixels). Notice that both algorithms leads to the same solution which gives a good reconstruction of

Denoising with Multiplicative noise

As imply by ( 8), the final estimate for each algorithm is given by taking the exponential of the result. In all experiments, Ψ was the ℓ1-norm. The Barbara image was set to a maximal intensity of 30 and the minimal to a non-zero value in order to avoid issues with the logarithm. The noise was added using M = 10 which leads to a medium level of noise. Fig 4 summarizes the results with the MAE and the execution times for Barbara, where the dictionary consisted of the curvelets transform. The difference between the visual results of the two algorithms can be explain be the fact we did not stop at the convergence, but instead stop after one thousand iterations. Our methods give correct reconstruction of the image. Fig. 5 displays the objective as a function of the iteration number and times (in s). Again, we can clearly see that Algorithm 2 converges faster than Algorithm 1. In this paper, we proposed two provably convergent algorithms for solving the linear inverse problems with a sparsity prior. The primal-dual proximal splitting algorithm seems to perform better in terms of convergence speed than the primal one. Moreover, its computational burden is lower than most comparable of state-of-art methods. Inverse problems with multiplicative noise does not enter currently in this framework, we will consider its adaptation to such problems in future work.

CONCLUSION

  is the class of all proper lower semicontinuous (lsc) convex functions from H to (-∞, +∞]. We denote by |||M||| = max x =0 Mx x the spectral norm of the linear operator M, and ker(M) := {x ∈ H : Mx = 0, x = 0} its kernel.

Algorithm 2 :

 2 Primal-dual scheme for solving (P f 1 ,γ,ψ ). Parameters: The observed image y, the dictionary Φ, number of iterations Niter, proximal steps σ > 0 and τ > 0, and regularization parameter γ > 0. Initialization: α0 = ᾱ0 = 0 ξ0 = η0 = 0. Main iteration: For t = 0 to Niter -1,• Data fidelity (Lemmas 8, 9 and 10):ξt+1 = (Iσ prox f 1 /σ )(ξt/σ + H • Φ ᾱt). • Positivity constraint: ηt+1 = (I -σPC)(ηt/σ + Φ ᾱt).• Sparsity-penalty (Lemma 11): αt+1 = prox τ γΨ `αtτ Φ T (H * ξt+1 + ηt+1)

Figure 1 :Figure 2 :

 12 Figure 1: Objective function for deconvolution under Poisson noise in function of iterations (left) and times (right).

Figure 3 :

 3 Figure 3: Objective function for inpainting with Gaussian noise in function of iterations (left) and times (right).

Figure 4 :

 4 Figure 4: Denoising results for Barbara using our two algorithms.

Figure 5 :

 5 Figure 5: Objective function for denoising with multiplicative noise in function of iterations (left) and times (right).

Table 1 :

 1 MAE and execution times for the deconvolution of the sky image.

		15]	RL-TV [10]	StabG [11]	PIDAL-FS [13]
	MAE	63.5	52.8	43	43.6
	Times	230s	4.3s	311s	342s
		Alg. 1	Alg. 2
	MAE	46	43.6
	Times	183s	154s


Diana Car
email: d.car@tue.nl

Jia Wang

Marcel A Verheijen

Erik P A M Bakkers

Sébastien R Plissard

Rationally Designed Single-Crystalline Nanowire Networks

Keywords: Nanowires, Nanocrosses, Networks, InSb, Growth Mechanisms

 allowing for configurations impossible to realize in thin film geometry.

Achieving a single crystalline nanowire junction is nevertheless a challenging task. [5,6] Defects are undesirable as they might induce carrier scattering and therefore reduce the transport properties. [7] We use the substrate crystallography to ensure the nanowires meet under the optimal angle for formation of single-crystalline structures.

Moreover, by exact control of the position of the catalyst in two directions we can tune the type, as well as the number, of wire-to-wire junctions.

Our method is generic and can be employed to synthesize interconnected nanowire architectures of group III-V, II-VI and IV materials. InSb is our material of choice due to its superior transport properties [8,9] which make it a promising platform for the future random access memories based on Majorana fermions (MFs). [10,11,12,13] The challenge ahead is to design and fabricate new quantum devices, which would allow for braiding of MFsan operation, which should bring us one step closer to topological quantum computing with these quasiparticles. The simplest circuit necessary to demonstrate the braiding of Majorana fermions consists of five superconducting islands connected by a 2-junction nanowire structure and enclosed in a transmission line resonator. [14] To fully exploit computational power of Majoranas, i.e. to allow for multiple independent braiding operations, complex nanowire systems with more than two junctions are demanded. [15] It is of great importance to develop a novel method of material synthesis, which would enable the controlled fabrication of planar interconnected nanowire architectures.

In order to control the nanowire position, arrays of gold islands are patterned on (001) InP substrate by electron beam lithography (EBL). InP nanowires grown on InP (001) substrates are used as stems for InSb nanowire growth to facilitate the nucleation of InSb wires, which is difficult on flat surfaces. [16,17] We use an InP substrate because P from the substrate cannot easily be incorporated into InSb due to the large miscibility gap in the ternary, [18] in contrast to As when an InAs substrate is used. Hence, a vertical <001> InP segment is first grown (Figure 1a). [19] It is then possible to change the growth direction of the InP nanowires into a <111>B direction (Figure 1b) by changing the contact angle of the catalyst droplet. [20] These kinked InP nanowires are then used as stems for the growth of InSb nanowires (Figure 1c). As already reported, [21,22,23] (001) substrates are ideal candidates for the growth of complex structures since two <111>B directions are available. To be more precise, on an (001) III-V substrate there are 4 <111> directions pointing upwards of which two correspond to <111>A and two to <111>B directions. As several studies have shown that InP and InSb nanowires normally grow in a <111>B direction, [24] we can conclude that there are only two growth directions available on an (001) III-V substrate (Figure 1d). Hence, 50% of the InSb wires will grow in the [11-1]B and 50% in the [-1-1-1]B direction. [20] Since the InSb wires are epitaxially connected to the substrate they will meet under an angle of 109.5° (Figure 1e), which is precisely the angle between two <111>B directions in a zincblende monocrystal. This, plus the fact that InSb wires crystallize in the defect-free zincblende crystal structure (as opposed to wurtzite InAs nanowires), is why our approach is expected to result in a high yield of single-crystalline junctions.

We first discuss the situation in which the catalyst particles are exactly aligned along the [-1-10] crystalline direction of the substrate (Figure 2a). When two InSb nanowires, growing towards each other (respectively in the [-1-1-1]B and the [11-1]B directions) meet, they touch tip-to-tip. The two gold droplets coalesce into a bigger one and a nanowire "bridge" is formed (Figure 2b,c). [21] Although this is an interesting approach to double the length of the InSb wires while maintaining the relatively small diameter, we are interested in multi-branched wires. In order to assist formation of crossing events, arrays with a misalignment in the y-direction were designed (Figure 2d). In this case, two InSb nanowires, growing towards each other from two mutually misaligned droplets, do not merge into a "bridge" structure. Instead, they grow in close proximity and finally merge due to radial nanowire growth forming a nanocross (Figure 2e,f). The crystalline structure of the junctions is investigated by transmission electron microscopy (TEM). Figure 3a shows a bright field TEM image of a crossed junction formed by two wires meeting under an optimal angle of 109.5°. The high resolution TEM images taken just below the droplets (Figure 3b,d) of two merged nanowires and the corresponding Fast-Fourier Transforms (FFT) (Figure 3c,e) show that the crystalline directions of both legs are identical. This structure is therefore a single crystal.

Nevertheless, the HR-TEM image of another junction (Figure S2, Supporting Information) shows a Moiré interference pattern which is a signature of two differently oriented crystals. The corresponding FFT pattern reveals a twinning diffraction pattern, with a 36.5° angle between the <200> directions of the two nanowires. This structure is not a single crystal, and there is a grain boundary at the junction. The formation of the grain boundary can be explained by rotational defects in the InP stem. The vertical <001> InP segment is defect-free whereas twinning occurs in the kinked <111>B InP segment (Figure S4, Supporting Information). After an odd number of rotational twins is formed in the <111>B InP nanowire, the top of the <111>B InP nanowire and the (001) InP substrate are twin-related. On the contrary, if there is an even number of rotational twins in the <111>B InP segment, its orientation is again identical to that of the substrate. Thus, if an even number of rotational twins is present in both the [11-1]B and [-1-1-1]B InP stem, the resulting crystalline orientations at the end of the two InP segments are identical. The same applies for InSb nanowires growing on top of them, and the resulting crossed junction is in that case single-crystalline.

It is important to notice that the orientation of a [11-1]B InP segment after an odd number of twins is not the same as the orientation of the [-1-1-1]B InP segment after an odd number of twins, since the rotations occur around different axes. In this case, a grain boundary is formed at the junction between the two InSb nanowires. The FFT pattern of such a twinned-twin junction will show a 36.5° angle between the two <200> directions, as shown in Figure S2, Supporting Information. Finally, an InSb nanocross can be formed from two InSb nanowires, one growing from an InP stem with an even number of twins, and one from an InP stem with an odd number of twins. The FFT diffraction pattern of such a twinned junction will show a 109.5° angle between the two <200> directions (Figure S3, Supporting Information). Taking into account all the possible crystalline orientations of InP nanowires, and the fact that there is a 50 % chance for an even and odd number of twin boundaries in an InP segment, we are left with four possible options: even-even, odd-odd, and even-odd/odd-even (Figure S5, Supporting Information); each of which has a 25% probability.

We note that with this method the yield of single crystalline junctions has been increased significantly compared with our previously reported approach (from 8% to 25%), [25] but further improvements are desirable for efficient device fabrication. One possible way to further increase the yield of single-crystalline nanostructures is to use defect-free InP nanowires as stems for InSb nanowire growth. In this manner we would be able to increase the yield of perfectly single-crystalline structures from 25% to 100%.

For this the growth of defect-free zinc-blende <111>B InP wires has to be developed, which has not been reported so far. An alternative route is to directly grow InSb nanowires on the vertical <001> InP nanowires, which are defect-free (Figure S4, Supporting Information). It is, however, difficult to nucleate InSb wires on top of these stems because during the cool down of the <001> InP nanowires a defect-free zincblende <111>B segment is formed. [20] The low surface energy of the AuInSb droplet on top of the kinked <111>B segment causes the wetting of the InP nanowire sidewalls, which will prevent InSb nanowires from nucleating.

A significant advantage of our method is the possibility of controlled fabrication of planar multiple wire structures. Important to mention is that the type of structure, which is predominantly formed, can be tuned by the catalyst spacing and offset. Figure 4a illustrates what we denote as a π-structure, as its shape resembles the one of the greek letter π. It is formed out of 3 wires, 2 of which are growing from mutually aligned catalyst particles and are therefore forming a nanowire-bridge structure, while the 3 rd nanowire is misaligned with respect to the other two and is crossing one of the two legs of a bridge. Figure 4b shows a slightly different πstructure, consisting of 3 wires and 2 crossed junctions. An example of a single-crystalline π -junction is shown in the Figure S1, Supporting Information, showing that this approach is suitable for fabricating key components of future Majorana random access memories. 4 nanowires can be merged in a way depicted in Figure 4c -two nano-bridges connected by one crossed junctionor they can form a closed nanowire loop as shown in Figure 4d. Bottom up synthesis of closed nanowire loops is an entirely novel concept.

A nanowire loop is formed out of four interconnected InSb nanowires, two of which are growing in the [11-1]B and the other two in [-1-1-1]B direction. In order to synthesize such complex structures one needs to be able to grow high yield <001> InP nanowires, exactly control the <001> to <111>B InP kinking event and be able to grow high aspect ratio InSb nanowires which should be at least 2.5 μm long in order to promote a double crossing event, necessary for the loop formation. Key concept here is again the misalignment of the gold droplets in an array, which allows us to grow crossed nanowire structures instead of nanowire-bridges. Our InSb nanowire loops are a promising material system, not only as building blocks of future Majorana random access memories, but also for exciting new physics experiments on coherent transport through the nanostructures from which the phase coherence length of electrons in InSb nanowires may be extracted.

Figure 4e illustrates a 6-wire structure: a closed 4 -nanowire loop whose one leg is forming a π-structure with the remaining 2 nanowires. A network formed out of 7 interconnected wires is shown in Figure 4f. Improvements in the length, as well as in the yield of InSb nanowires on (001) substrate, would allow for controllable growth of even larger InSb nanowire networks. Furthermore, the method presented is general, such that it can be applied to assemble complex nanowire structures of other group III-V, II-VI and IV semiconductor materials and even hetero materials.

To summarize, we demonstrate an unprecedented level of control in growth of complex InSb nanowire structures. We report, for the first time, planar, crystalline nanowire structures formed out of more than two nanowires -nanowire π-structures, closed nanowire loops and nanowire networks The perfect control of the droplet positioning, i.e. the offset between the neighboring gold droplets in an array, allows us to favor crossed junctions over "bridge"-junctions and, in such a way, increase the yield of advantageous structures.

Experimental Section

Nanowire growth: Nanowires have been grown on (001) InP substrate using Aucatalyzed vapor-liquid-solid process by metal organic vapor phase epitaxy in an Aixtron MOVPE machine. InP <001> nanowires were first grown at a temperature of 440 °C for 9 min using tri-methyl-indium (TMI) and phosphine (PH 3 ) as precursors with growth pressure of 50 mbar and precursor molar fractions Xi (TMI) = 2.42 x 10 -5 and Xi (PH 3 ) = 7.44 x 10 -3 . To increase the yield of vertical wires growing in the <001> direction the Au droplet shape and composition was optimized prior to wire growth: when the system reached the temperature of 440 °C under PH 3 pressure, we simultaneously switched PH 3 off and TMI on for 15 s. These high yield <001> InP nanowires are used as templates for switching of the growth direction. After 9 min of <001> nanowire growth, TMI is switched off for 2 min while PH 3 is kept on. When we re-start the growth by switching TMI on again, nanowires grow in the <111>B direction. InP <111>B segments are grown using the same temperature and precursor flows as the <001> InP wires.

These kinked InP nanowires are then taken out of the growth reactor and transported, in nitrogen atmosphere, to the other Aixtron MOVPE machine dedicated to the growth of antimonides only. InSb nanowires are grown on top of the kinked InP stems at 495 °C In 25% of the cases, the nanostructure will be single crystalline. A defect-free pure zinc blende (111) InP segment would lead to a 100% yield of single crystalline junctions.

Figure 1 .

 1 Figure 1. InSb nanowires grown on an (001) InP substrate. (a) <001> InP nanowires

Figure 2 .

 2 Figure 2. Positionand directioncontrolled growth of InSb nanowire networks.

Figure 3 .

 3 Figure 3. Transmission Electron Microscopy (TEM) study of a single-crystalline

Figure 4 .

 4 Figure 4. Nanowire networks. (a) 30° tilted SEM image of the nanowire πstructure.

Figure S1 :

 S1 Figure S1: Single-crystalline InSb πstructure.

Figure S2 :

 S2 Figure S2: InSb nanocross with a twinned-twin junction.

  (d) FFT corresponding to HRTEM image in (b). There is a twinned-twin relation between crystallographic orientations of the two nanowires, the angle between the two <200> directions being 36.5°.

Figure S3 :

 S3 Figure S3: InSb T-structure with a twinned junction.

Table S5 :

 S5 Table summarizing the four possible crossing scenarios.

for 20 min, using tri-methyl-indium (TMI) (Xi= 1.1 x 10 -5 ) and tri-methyl-antimony (TMSb) (Xi= 3.4 x 10 -3 ).

Acknowledgements

This work has been supported by the Dutch Organization for Scientific Research (NWO), the Foundation for Fundamental Research on Matter (FOM) and the European Union Seventh Framework Programme under grant agreement no. 265073.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.

Additional information

The authors declare no competing financial interests.

Supporting Information

Rationally designed single-crystalline nanowire networks

Diana Car 1,* , Jia Wang 1 , Marcel A. Verheijen 1,2 , Erik P.A.M. Bakkers 1,3 , Sébastien R. Plissard 3,4,5 Summary of the supplementary materials: From the bright field TEM image we can see that the <001> InP segment and <111> InSb nanowire are defect free. On the contrary, a significant number of stacking faults present in the <111> InP segment makes it safe to conclude that there is a 50:50 chance an odd or even number of twin domains will be present in this segment.