Elvis Dohmatob

Arthur Mensch

Gael Varoquaux

Bertrand Thirion

Learning brain regions via large-scale online structured sparse dictionary-learning

) reference method, the online nature of our proposed algorithm allows it to scale to arbitrarily sized datasets. Experiments on brain data show that our proposed method extracts structured and denoised dictionaries that are more intepretable and better capture inter-subject variability in small medium, and large-scale regimes alike, compared to state-of-the-art models.

Introduction

In neuro-imaging, inter-subject variability is often handled as a statistical residual and discarded. Yet there is evidence that it displays structure and contains important information. Univariate models are ineffective both computationally and statistically due to the large number of voxels compared to the number of subjects. Likewise, statistical analysis of weak effects on medical images often relies on defining regions of interests (ROIs). For instance, pharmacology with Positron Emission Tomography (PET) often studies metabolic processes in specific organ sub-parts that are defined from anatomy. Population-level tests of tissue properties, such as diffusion, or simply their density, are performed on ROIs adapted to the spatial impact of the pathology of interest. In functional brain imaging, e.g function magnetic resonance imaging (fMRI), ROIs must be adapted to the cognitive process under study, and are often defined by the very activation elicited by a closely related process [START_REF] Saxe | Divide and conquer: a defense of functional localizers[END_REF]. ROIs can boost statistical power by reducing multiple comparisons that plague image-based statistical testing. If they are defined to match spatially the differences to detect, they can also improve the signal-to-noise ratio by averaging related signals. However, the crux of the problem is how to define these ROIs in a principled way. Indeed, standard approaches to region definition imply a segmentation step. Segmenting structures on first-level statistical maps, as in fMRI, typically yields meaningful units, but is limited by the noise inherent to these maps. Relying on a different imaging modality hits cross-modality correspondence problems.

Sketch of our contributions.

In this manuscript, we propose to use the variability of the statistical maps across the population to define regions. This idea is reminiscent of clustering approaches, that have been employed to define spatial units for quantitative analysis of information as diverse as brain fiber tracking, brain activity, brain structure, or even imaging-genetics. See [START_REF] Varol | Supervised block sparse dictionary learning for simultaneous clustering and classification in computational anatomy[END_REF][START_REF] Hibar | Genetic clustering on the hippocampal surface for genome-wide association studies[END_REF] and references therein. The key idea is to group together features -voxels of an image, vertices on a mesh, fiber tracts-29th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

based on the quantity of interest, to create regions -or fiber bundles-for statistical analysis. However, unlike clustering that models each observation as an instance of a cluster, we use a model closer to the signal, where each observation is a linear mixture of several signals. The model is closer to mode finding, as in a principal component analysis (PCA), or an independent component analysis (ICA), often used in brain imaging to extract functional units [START_REF] Beckmann | Probabilistic independent component analysis for functional magnetic resonance imaging[END_REF]. Yet, an important constraint is that the modes should be sparse and spatially-localized. For this purpose, the problem can be reformulated in a dictionary-learning, as ICA/PCA, with appropriate spatial and sparse penalties [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF][START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF].

We propose a multivariate online dictionary-learning method for obtaining decompositions with structured and sparse components (aka atoms). Sparsity is to be understood in the usual sense: the atoms contain mostly zeros. This is imposed via an 1 penalty on the atoms. By "structured", we mean that the atoms are piece-wise smooth and compact, thus making up blobs, as opposed to scattered patterns of activation. We impose this type of structure via a Laplacian penalty on the dictionary atoms. Combining the two penalties, we therefore obtain decompositions that are closer to known functional organization of the brain. This non-trivially extends the online dictionary-learning / dictionary-learning work [START_REF]Online learning for matrix factorization and sparse coding[END_REF], with only a factor of 2 or 3 on the running time. By means of experiments on a large public dataset, we show the improvements brought by the spatial regularization with respect to traditional 1 -regularized dictionary learning. We also provide a concise study of the impact of hyper-parameter selection on this problem and describe the optimality regime, based on relevant criteria (reproducibility, captured variability, explanatory power in prediction problems).

Smooth Sparse Online Dictionary-Learning (Smooth-SODL)

Consider a stack X ∈ R n×p of n subject-level brain images X 1 , X 2 , . . . , X n each of shape n 1 × n 2 × n 3 , seen as p-dimensional row vectors -with p = n 1 × n 2 × n 3 , the number of voxels. These could be images of fMRI activity patterns like statistical parametric maps of brain activation, raw pre-registered (into a common coordinate space) fMRI time-series, PET images, etc. We would like to decompose these images as a product of k ≤ min(n, p) component maps (aka latent factors or dictionary atoms)

V 1 , . . . , V k ∈ R p×1 and modulation coefficients U 1 , . . . , U n ∈ R k×1 called codes (one k-dimensional code per sample point), i.e X i ≈ VU i , for i = 1, 2, . . . , n (1)
where

V := [V 1 | . . . |V k ] ∈ R p×k
, an unknown dictionary to be estimated. Typically, p ∼ 10 5 -10 6 (in full-brain high-resolution fMRI) and n ∼ 10 2 -10 5 (for example, in considering all the 500 subjects and all the about 15 -20 functional tasks of the Human Connectome Project dataset [START_REF] Essen | The Human Connectome Project: A data acquisition perspective[END_REF]). Our work handles the extreme case where both n and p are large (massive-data setting). It is reasonable then to only consider under-complete dictionaries: k ≤ min(n, p). Typically, we use k ∼ 50 or 100 components. It should be noted that online optimization is not only crucial in the case where n/p is big; it is relevant whenever n is large, leading to prohibitive memory issues irrespective of how big or small p is.

As explained in the introduction, we would want the component maps (aka dictionary atoms) V j to be sparse and spatially smooth. A principled way to achieve such a goal is to impose a boundedness constraint on 1 -like norms of these maps to achieve sparsity and simultaneously impose smoothness by penalizing their Laplacian. Thus, we propose the following penalized dictionary-learning model

min V∈R p×k lim n→∞ 1 n n i=1 min Ui∈R k 1 2 X i -VU i 2 2 + 1 2 α U i 2 2 + γ k j=1 Ω Lap (V j ). subject to V 1 , . . . , V k ∈ C (2)
The ingredients in the model can be broken down as follows:

• Each of the terms max Ui∈R k 1 2 X i -VU i 2
2 measures how well the current dictionary V explains data X i from subject i. The Ridge penalty term φ(U i ) ≡ 1 2 α U i 2 2 on the codes amounts to assuming that the energy of the decomposition is spread across the different samples. In the context of a specific neuro-imaging problem, if there are good grounds to assume that each sample / subject should be sparsely encoded across only a few atoms of the dictionary, then we can use the 1 penalty φ(U i ) := α U i 1 as in [START_REF]Online learning for matrix factorization and sparse coding[END_REF]. We note that in contrast to the 1 penalty, the Ridge leads to stable codes. The parameter α > 0 controls the amount of penalization on the codes.

• The constraint set C is a sparsity-inducing compact simple 1

Ω Lap (v) := 1 2 p k=1 (∇ x v) 2 k + (∇ y v) 2 k + (∇ z v) 2 k = 1 2 v T ∆v ≥ 0, ∀v ∈ R p , (3)
∇ x being the discrete spatial gradient operator along the x-axis (a p-by-p matrix), ∇ y along the y-axis, etc., and ∆ := ∇ T ∇ is the p-by-p matrix representing the discrete Laplacian operator. This penalty is meant to impose blobs. The regularization parameter γ ≥ 0 controls how much regularization we impose on the atoms, compared to the reconstruction error.

The above formulation, which we dub Smooth Sparse Online Dictionary-Learning (Smooth-SODL) is inspired by, and generalizes the standard dictionary-learning / dictionary-learning framework of [START_REF]Online learning for matrix factorization and sparse coding[END_REF] -henceforth referred to as Sparse Online Dictionary-Learning (SODL); setting γ = 0, we recover SODL [START_REF]Online learning for matrix factorization and sparse coding[END_REF].

3 Estimating the model

Algorithms

The objective function in problem ( 2) is separately convex and block-separable w.r.t each of U and V but is not jointly convex in (U, V). Also, it is continuously differentiable on the constraint set, which is compact and convex. Thus by classical results (e.g Bertsekas [START_REF] Bertsekas | Nonlinear programming[END_REF]), the problem can be solved via Block-Coordinate Descent (BCD) [START_REF]Online learning for matrix factorization and sparse coding[END_REF]. Reasoning along the lines of [START_REF] Jenatton | Structured sparse principal component analysis[END_REF], we derive that the BCD iterates are as given in Alg. 1. A crucial advantage of using a BCD scheme is that it is parameter free: there is not step size to tune. The resulting algorithm Alg. 1, is adapted from [START_REF]Online learning for matrix factorization and sparse coding[END_REF]. It relies on Alg. 2 for performing the structured dictionary updates, the details of which are discussed below.

Algorithm 1 Online algorithm for the dictionary-learning problem (2)

Require: Regularization parameters α, γ > 0; initial dictionary V ∈ R p×k , number of passes / iterations T on the data.

1: A 0 ← 0 ∈ R k×k , B 0 ← 0 ∈ R p×k (historical "sufficient statistics") 2: for t = 1 to T do 3:
Empirically draw a sample point X t at random.

4:

Code update: Ridge-regression (via SVD of current dictionary V)

U t ← argmin u∈R k 1 2 X t -Vu 2 2 + 1 2 α u 2 2 . (4) 5:
Rank-1 updates:

A t ← A t-1 + U t U T t , B t ← B t-1 + X t U T t 6:
BCD dictionary update: Compute update for dictionary V using Alg. 2. 7: end for Update of the codes: Ridge-coding. The Ridge sub-problem for updating the codes

U t = (V T V + αI) -1 V T X t (5)
is computed via an SVD of the current dictionary V. For α ≈ 0, U t reduces to the orthogonal projection of X t onto the image of the current dictionary V. As in [START_REF]Online learning for matrix factorization and sparse coding[END_REF], we speed up the overall algorithm by sampling mini-batches of η samples X t , . . . , X η and compute the corresponding codes U 1 , U 2 , ..., U η at once. We typically use we use mini-batches of size η ∼ 20.

1 Mainly in the sense that the Euclidean projection onto C should be easy to compute.

3

BCD dictionary update for the dictionary atoms. Let us define time-varying matrices

A t := t i=1 U i U T i ∈ R k×k and B t := t i=1 X i U T i ∈ R p×k
, where t = 1, 2, . . . denotes time. We fix the matrix of codes U, and for each j, consider the update of the jth dictionary atom, with all the other atoms V k =j kept fixed. The update for the atom V j can then be written as

V j = argmin v∈C,V=[V 1 |...|v|...|V k ] 1 t t i=1 1 2 X i -VU i 2 2 + γΩ Lap (v) = argmin v∈C,V=[V 1 |...|v|...|V k ] t i=1 1 2 X i -VU i 2 2 + γtΩ Lap (v) = argmin v∈C F γ(At[j,j]/t) -1 (v, V j + A t [j, j] -1 (B j t -VA j t )
refer to [START_REF]Online learning for matrix factorization and sparse coding[END_REF] for the details

), (6)
where

F γ (v, a) ≡ 1 2 v -a 2 2 + γΩ Lap (v) = 1 2 v -a 2 2 + 1 2 γv T ∆v.
Algorithm 2 BCD dictionary update with Laplacian prior

Require: V = [V 1 | . . . |V k ] ∈ R p×k (input dictionary), 1: A = [A 1 | . . . |A k ] ∈ R k×k , B t = [B 1 t | . . . |B k t ] ∈ R p×k (
history) 2: while stopping criteria not met, do 3:

for j = 1 to r do 4:
Fix the code U and all atoms k = j of the dictionary V and then update V j as follows

V j ← argmin v∈C F γ(At[j,j]/t) -1 (v, V j + A t [j, j] -1 (B j t -VA j )) (7)
(See below for details on the derivation and the resolution of this problem) 5:

end for 6: end while Problem ( 6) is thus a compactly-constrained minimization of the 1-strongly-convex quadratic functions F γ (., a) : R p → R defined above. This problem can further be identified with a denoising instance (i.e in which the design matrix or deconvolution operator is the identity operator) of the GraphNet model [START_REF] Grosenick | Interpretable whole-brain prediction analysis with GraphNet[END_REF][START_REF] Hebiri | The Smooth-Lasso and other 1 + 2 -penalized methods[END_REF]. Fast first-order methods like FISTA [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF] with optimal rates O(L/ √ ) are available2 for solving such problems to arbitrary precision > 0. One computes the Lipschitz constant to be L Fγ (.,a) ≡ 1 + γL ΩLap = 1 + 4Dγ, where as before, D is the number of spatial dimensions with D = 3 for volumic images. One should also mention that under certain circumstances, it is possible to perform the dictionary updates in the Fourier domain, via FFT. This alternative approach is developed in Appendix B.

Finally, one notes that, since constraints in problem (2) are separable in the dictionary atoms V j , the BCD dictionary-update algorithm Alg. 2 is guaranteed to converge to a global optimum, at each iteration [START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF]Online learning for matrix factorization and sparse coding[END_REF].

How difficult is the dictionary update for our proposed model ? A favorable property of the vanilla dictionary-learning [START_REF]Online learning for matrix factorization and sparse coding[END_REF] is that the BCD dictionary updates amount to Euclidean projections onto the constraint set C, which can be easily computed for a variety of choices (simplexes, closed convex balls, etc.). One may then ask: do we retain a comparable algorithmic simplicity even with the additional Laplacian terms Ω Lap (V j ) ? The short answer is yes: empirically, we found that 1 or 2 iterations of FISTA [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF] are sufficient reach an accuracy of 10 -6 in problem [START_REF] Bertsekas | Nonlinear programming[END_REF], which is sufficient to obtain a good decomposition in the overall algorithm.

However, choosing γ "too large" will provably cause the dictionary updates to eventually take forever to run. Indeed, the Lipschitz constant in problem ( 6) is L t = 1 + 4Dγ(A t [j, j]/t) -1 , which will blow-up (leading to arbitrarily small step-sizes) unless γ is chosen so that

γ = γ t = O max 1≤j≤k A t [j, j] = O max 1≤j≤k t i=1 U j 2 2 /t = O( A t ∞,∞ /t). (8)
Finally, the Euclidean projections onto the 1 ball C can be computed exactly in linear-time O(p) (see for example [START_REF] Condat | Fast projection onto the simplex and the 1 ball[END_REF][START_REF] Duchi | Efficient projections onto the l 1-ball for learning in high dimensions[END_REF]). The dictionary atoms j are repeatedly cycled and problem (6) solved. All in all, in practice we observe that a single iteration is sufficient for the dictionary update sub-routine in Alg.

2 to converge to a qualitatively good dictionary. Hyper-parameter tuning. Parameterselection in dictionary-learning is known to be a difficult unsolved problem [START_REF]Online learning for matrix factorization and sparse coding[END_REF][START_REF] Jenatton | Structured sparse principal component analysis[END_REF], and our proposed model ( 2) is not an exception to this rule. We did an extensive study of the quality of estimated dictionary varies with the model hyper-parameters (α, γ, τ ). The data experimental setup is described in Section 5. The results are presented in Fig. 1. We make the following observations: Taking the sparsity parameter τ in (2) too large leads to dense atoms that perfectly explain the data but are not very intepretable. Taking it too small leads to overly sparse maps that barely explain the data. This normalized sparsity metric (small is better, ceteris paribus) is defined as the mean ratio V j 1 / V j 2 over the dictionary atoms. Concerning the α parameter, inspired by [START_REF] Ying | Online regularized classification algorithms[END_REF], we have found the following time-varying data-adaptive choice for the α parameter to work very well in practice:

Convergence

α = α t ∼ t -1/2 . ( 9
)
Likewise, care must be taken in selecting the Laplacian regularization parameter γ. Indeed taking it too small amounts to doing vanilla dictionary-learning model [START_REF]Online learning for matrix factorization and sparse coding[END_REF]. Taking it too large can lead to degenerate maps, as the spatial regularization then dominates the reconstruction error (data fidelity) term. We find that there is a safe range of the parameter pair (γ, τ ) in which a good compromise between the sparsity of the dictionary (thus its intepretability) and its explanation power of the data can be reached. See Fig. 1. K-fold cross-validation with explained variance metric was retained as a good strategy for setting the Laplacian regularization γ parameter and the sparsity parameter τ .

Initialization of the dictionary. Problem ( 2) is non-convex jointly in (U, V), and so initialization might be a might be a crucial issue. However, in our experiments, we have observed that even randomly initialized dictionaries eventually produce sensible results that do not jitter much across different runs of the same experiment.

Related works

While there exist algorithms for online sparse dictionary-learning that are very efficient in large-scale settings (for example [START_REF]Online learning for matrix factorization and sparse coding[END_REF], or more recently [START_REF] Mensch | Dictionary Learning for Massive Matrix Factorization[END_REF]) imposing spatial structure introduces couplings in the corresponding optimization problem [START_REF] Dohmatob | Benchmarking solvers for TV-l1 least-squares and logistic regression in brain imaging[END_REF]. So far, spatially-structured decompositions have been solved by very slow alternated optimization [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF][START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF]. Notably, structured priors such as TV-L1

[3] minimization, were used by [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF] to extract data-driven state-of-the-art atlases of brain function. However, alternated minimization is very slow, and large-scale medical imaging has shifted to online solvers for dictionary-learning like [START_REF]Online learning for matrix factorization and sparse coding[END_REF] and [START_REF] Mensch | Dictionary Learning for Massive Matrix Factorization[END_REF]. These do not readily integrate structured penalties. As a result, the use of structured decompositions has been limited so far, mostly due to the computational cost of the ensuing algorithms. Our approach instead uses a Laplacian penalty to impose spatial structure at a very minor cost and adapts the online-learning dictionary-learning framework [START_REF]Online learning for matrix factorization and sparse coding[END_REF], resulting in a fast and scalable structured decomposition. Second, the approach in [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF] though very novel, is heuristic, as it does not come with theoretical guarantees. In contrast, our method enjoys the same convergence guarantees and comparable numerical complexity as the basic unstructured online dictionary-learning [START_REF]Online learning for matrix factorization and sparse coding[END_REF].

Finally, one should also mention [START_REF] Varoquaux | Cohort-level brain mapping: learning cognitive atoms to single out specialized regions[END_REF] which introduced an online group-level functional brain mapping strategy for differentiating regions reflecting the variety of brain network configurations observed a the population, by learning a sparse-representation of these in the spirit of [START_REF]Online learning for matrix factorization and sparse coding[END_REF].

Experiments

Setup. Our experiments were done on task fMRI data from 500 subjects from the HCP -Human Connectome Project-dataset [START_REF] Essen | The Human Connectome Project: A data acquisition perspective[END_REF]. These task fMRI data were acquired in an attempt to assess major domains that are thought to sample the diversity of neural systems of interest in functional connectomics. We studied the activation maps related to a task that involves language (story understanding) and mathematics (mental computation). This particular task is expected to outline number, attentional and language networks, but the variability modes observed in the population cover even wider systems.

For the experiments, mass-univariate General Linear Models (GLMs) [START_REF] Friston | Statistical Parametric Maps in Functional Imaging: A General Linear Approach[END_REF] for n = 500 subjects were estimated for the Math vs Story contrast (language protocol), and the corresponding full-brain Z-score maps each containing p = 2.6 × 10 5 voxels, were used as the input data X ∈ R n×p , and we sought a decomposition into a dictionary of k = 40 atoms (components). The input data X were shuffled and then split into two groups of the same size.

Models compared and metrics. We compared our proposed Smooth-SODL model ( 2) against both the Canonical ICA -CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF], a single-batch multi-subject PCA/ICA-based method, and the standard SODL (sparse online dictionary-learning) [START_REF]Online learning for matrix factorization and sparse coding[END_REF]. While the CanICA model accounts for subject-to-subject differences, one of its major limitations is that it does not model spatial variability across subjects. Thus we estimated the CanICA components on smoothed data: isotropic FWHM of 6mm, a necessary preprocessing step for such methods. In contrast, we did no pre-smoothing for the SODL of Smooth-SODL models. The different models were compared across a variety of qualitative and quantitative metrics: visual quality of the dictionaries obtained, explained variance, stability of the dictionary atoms, their reproducibility, performance of the dictionaries in predicting behavioral scores (IQ, picture vocabulary, reading proficiency, etc.) shipped with the HCP data [START_REF] Essen | The Human Connectome Project: A data acquisition perspective[END_REF]. For both SODL [START_REF]Online learning for matrix factorization and sparse coding[END_REF] and our proposed Smooth-SODL model, the constraint set for the dictionary atoms was taken to be a simplex C := S p (τ ) (see section 2 for definition). The results of these experiments are presented in Fig. 2 and Tab. 1.

Results

Running time. On the computational side, The standard dictionary-learning SODL algorithm [START_REF]Online learning for matrix factorization and sparse coding[END_REF] with a batch size of η = 20 took about 110s (≈ 1.7 minutes) to run, whilst with the same batch size, our proposed Smooth-SODL model ( 2) implemented in Alg. 1 took 340s (≈ 5.6 minutes), which is slightly less than 3× slower than SODL. Finally, CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF] for this experiment took 530s (≈ 8.8 minutes) to run, which is about 5× slower than the SODL model and 1.6× slower than our proposed Smooth-SODL (2) model. All experiments were run on a single CPU of a modern laptop.

Qualitative assessment of dictionaries. As can be seen in Fig. 2(a), all methods recover dictionary atoms that represent known functional brain organization; notably the dictionaries all contain the well-known executive control and attention networks, at least in part. Vanilla dictionary-learning leverages the denoising properties of the 1 sparsity constraint, but the voxel clusters are not very structured. For, example most blobs are surrounded with a thick ring of very small nonzero values. In contrast, our proposed regularization model leverages both sparse and structured dictionary atoms, which are more spatially structured and less noisy.

In contrast to both SODL and Smooth-SODL, CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF] is an ICA-based method which enforces no notion of sparsity whatsoever. The result are therefore dense and noisy dictionary atoms that explain the data very well (Fig. 2(b) but which are completely unintepretable. In a futile attempt to remedy the situation, in practice such PCA/ICA-based methods (including FSL's MELODIC tool (a) Qualitative comparison of the estimated dictionaries. Each column represents an atom of the estimated dictionary, where atoms from the different models (the rows of the plots) have been matched via a Hungarian algorithm. Here, we only show a limited number of the most "intepretable" atoms. Notice how the major structures in each atom are reproducible across the different models. Maps corresponding to hard-thresholded CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF] components have also been included, and have been called tCanICA. In contrast, the maps from the SODL [START_REF]Online learning for matrix factorization and sparse coding[END_REF] and our proposed Smooth-SODL (2) were not been thresholded. 2) model against competing state-of-the-art methods like SODL (sparse online dictionary-learning) [START_REF]Online learning for matrix factorization and sparse coding[END_REF] and CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF].

[21]) are hard-thresholded in order to see information. For CanICA, the hard-thresholded version has been named tCanICA in Fig. 2. That notwithstanding, notice how the major structures (parietal lobes, sulci, etc.) in each atom are reproducible across the different models.

Stability-fidelity trade-offs. PCA/ICA-based methods like CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF] and MELODIC [START_REF] Smith | Advances in functional and structural MR image analysis and implementation as FSL[END_REF] are the optimal linear decomposition method to maximize explained variance on a dataset. On the training set, CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF] outperforms all others algorithms with about 66% (resp. 50% for SODL [START_REF]Online learning for matrix factorization and sparse coding[END_REF] and 58% for Smooth-SODL) of explained variance on the training set, and 60% (resp. 49% for SODL and 55% for Smooth-SODL) on left-out (test) data. See Fig. 2(b). However, as noted in the above paragraph, such methods lead to dictionaries that are hardly intepretable and thus the user must recourse to some kind of post-processing hard-thresholding step, which destroys the estimated model. More so, assessing the stability of the dictionaries, measured by mean correlation between corresponding atoms, across different splits of the data, CanICA [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF] scores a meager 0.1, whilst the hard-thresholded version tCanICA obtains 0.2, compared to 0.4 for Smooth-SODL and 0.1 for SODL.

Is spatial regularization really needed ?

As rightly pointed out by one of the reviewers, one does not need spatial regularization if data are abundant (like in the HCP). So we computed learning curves of mean explained variance (EV) on test data, as a function of the amount training data seen by both Smooth-SODL and SODL [START_REF]Online learning for matrix factorization and sparse coding[END_REF] (Table 1). In the beginning of the curve, our proposed spatially regularized Smooth-SODL model starts off with more than 31% explained variance (computed on 241 subjects), after having pooled only 17 subjects. In contrast, the vanilla SODL model [START_REF]Online learning for matrix factorization and sparse coding[END_REF] scores a meager 2% explained variance; this corresponds to a 14-fold gain of Smooth-SODL over SODL. As more and more that are pooled, both models explain more variance, and the gap between Smooth-SODL and SODL reduces, and both models perform comparably asymptotically.

Nb. subjects pooled mean EV for vanilla SODL [START_REF]Online learning for matrix factorization and sparse coding[END_REF] Thus our proposed Smooth-SODL method extracts structured denoised dictionaries that better capture inter-subject variability in small, medium, and large-scale regimes alike.

Prediction of behavioral variables. If Smooth-SODL captures the patterns of inter-subject variability, then it should be possible to predict cognitive scores y like picture vocabulary, reading proficiency, math aptitude, etc. (the behavioral variables are explained in the HCP wiki [START_REF]HCP wiki[END_REF]) by projecting new subjects' data into this learned low-dimensional space (via solving the ridge problem (4) for each sample X t ), without loss of performance compared with using the raw Z-values values X. Let RAW refer to the direct prediction of targets y from X, using the top 2000 most voxels most correlated with the target variable. Results of for the comparison are shown in Fig. 2(c). Only variables predicted with a a positive mean (across the different methods and across subjects) R-score are reported. We see that the RAW model, as expected over-fits drastically, scoring an R 2 of 0.3 on training data and only 0.14 on test data. Overall, for this metric CanICA performs best than all the other models in predicting the different behavioral variables on test data. However, our proposed Smooth-SODL model outperforms both SODL [START_REF]Online learning for matrix factorization and sparse coding[END_REF] and tCanICA, the thresholded version of CanICA.

Concluding remarks

To extract structured functionally discriminating patterns from massive brain data (i.e data-driven atlases), we have extended the online dictionary-learning framework first developed in [START_REF]Online learning for matrix factorization and sparse coding[END_REF], to learn structured regions representative of brain organization. To this end, we have successfully augmented [START_REF]Online learning for matrix factorization and sparse coding[END_REF] with a Laplacian prior on the component maps, while conserving the low numerical complexity of the latter. Through experiments, we have shown that the resultant model -Smooth-SODL model ( 2)extracts structured and denoised dictionaries that are more intepretable and better capture inter-subject variability in small medium, and large-scale regimes alike, compared to state-of-the-art models. We believe such online multivariate online methods shall become the de facto way do dimensionality reduction and ROI extraction in future.

Implementation. The authors' implementation of the proposed SS0MF (2) model will soon be made available as part of the Nilearn package [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF].

A Convergence of the proposed algorithm

We now show how the convergence of our proposed algorithm follows effortless from [START_REF]Online learning for matrix factorization and sparse coding[END_REF]. Note the the objective in (2) can be rewritten as function of the dictionary alone like so

E x ( (x, V)) + γ k j=1 Ω Lap (V j ) a.s = lim t→∞ f t (V) + γ k j=1 Ω Lap (V j ),
where (x, V) := min u∈R k 1 2 x -Vu 2 2 + 1 2 α u 2 2 , and f t (V) := 1 t t i=1 (X i , V), with the X i 's sampled from the data. For each time t ≥ 0, define

ft (V) := 1 t t i=1 1 2 X i -V Ûi 2 2 + 1 2 α Ûi 2 2 + γ k j=1 Ω Lap (V j ) (10)
with each code Ûi is computed online by solving the Ridge problem (4). The following observations are immediate:

• Assumption (A) of [START_REF]Online learning for matrix factorization and sparse coding[END_REF], that the data distribution admits compact support, is automatically true for all MRI and PET data, because it is imposed by acquisition (for example, there is no data beyond the bounding-box).

• Assumption (C), that the solution to the coding problem, which corresponds to (4) in our case, is unique and Lipschitz continuous w.r.t to the incoming data X t is automatically satisfied for us since Ût (X t ) = (V T V + αI) -1 V T X t , a linear transformation of X t .

As a consequence, we have the classical convergence guarantees as in in [START_REF]Online learning for matrix factorization and sparse coding[END_REF]:

• Almost-sure convergence of the dictionary: V(t) -V(t -1) = O(1/t) a.s. as t → ∞.

• Almost-sure convergence of the risk:

ft (V(t)) converges a.s. as t → ∞.

-f t (V(t)) -f (V(t)) → 0 a.s. as t → ∞.

-f t (V(t)) converges a.s. as t → ∞.

Finally, one notes that since constraints in problem (2) are separable in the dictionary atoms V j , the BCD dictionary-update algorithm Alg. 2 is guaranteed to converge to a global optimum. [START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF]Online learning for matrix factorization and sparse coding[END_REF].

B Working in frequency domain, when it is possible.

To close this section, let us point out a few special instances cases of problem [START_REF] Bertsekas | Nonlinear programming[END_REF], for peculiar choices of the constraint set Q. First note that the objective in problem ( 6) can be conveniently rewritten as

F γAt[j,j] -1 (v, V j + A t [j, j] -1 (VA j -B j t )) = 1 2 (v -Ṽj ) T (I -γA t [j, j] -1 ∆)(v -Ṽj ) = 1 2 (v -Vj ) T (I -γA t [j, j] -1 ∆)(v -Vj ), (11)
with

Ṽj := (A t [j, j]I -γ∆) -1 V j + A t [j, j] -1 (VA j -B j t ) . (12)
We note that the matrix-inversion (I -γ∆) -1 that appears in the formula above is a Laplacian filter, and can be efficiently applied in closed-form (i.e non-iteratively) in the Fourier / frequency domain. Indeed, under periodic boundary conditions, the discrete Laplacian ∆ is Block-Circulant with Circulant Blocks (BCCB) and so is diagonalizable in the Fourier domain. Precisely,

∆ = F * ΛF (13)
where the complex orthonormal operator F represents the fast Discrete Fourier Transform (DFT), and Λ is diagonal matrix made p eigenvalues (including multiplicities) of the Laplace operator ∆, given by

Λ(ω) := - 3 d=1 2 sin ω d π 2n d 2 = -2 3 d=1 1 -cos ω d π n d ≤ 0, for ω = (ω 1 , ω 2 , ω 3 ) ∈ [[0, n 1 -1]] × [[0, n 2 -1]] × [[0, n 3 -1]].
We note that the spectral norm of Laplace operator in D dimensions (here

D = 3) is ∆ 2 = γmax (-∆) = 2 × D × (1 + 1) = 4D.
Now, one can then harvest the closed-form solution

(I -γ∆) -1 a = (F -1 (I -γΛ) -1 F)(a) = F -1 (s) , (14)
where s ∈ R p is defined by s(ω) := â(ω)

1 -γ ∆(ω) , with â := F(a). These DFT computations have complexity O(p log p).

For applying the DFTs above, one can use the FFTW library for computing the forward and inverse Fourier transforms needed to apply the Laplacian filter (refer to paragraph B). FFTW is generally taught to be one of the fastest implementations of the FFT, yielding up to 3× speedup against competing libraries like LAPACK.

Pure 2 constraint. Here, the constraint set C is an L2 ball (with radius = 1, w.l.o.g) in R 2 . By the Rayleigh energy theorem (aka Parseval's identity for the DFT), one has

v 2 = p v 2 2
, ∀v ∈ R p and so problem (6) can be written as

V j ← argmin v∈R p , v 2
where

E := (I -γA t [j, j] -1 Λ) 1 2 v s.t v ∈ C p , v 2 2 ≤ p , (16)
a hyper-ellipsoid in standard position (i.e 0-centered and axes-aligned). Using elementary geometric arguments, one can show that the projection P E ( Vj ) can be computed efficiently using a kind of root-finding algorithm [START_REF] Dai | Fast Algorithms for Projection on an Ellipsoid[END_REF], and converges exponentially fast.

Non-negative Lasso. In case the constraint set C for the dictionary atoms is a simplex S p (τ ), the simplex (see section 2), then the BCD update for the jth atom becomes

V j ← argmin v∈R p , v≥0, 1 T v≤1 1 2 (v -Vj ) * (I -γA t [j, j] -1 Λ)(v -Vj ) = F * argmin v∈C p , -F * v≤0, 1T v≤1 1 2 (v -Vj ) * (I -γA t [j, j] -1 Λ)(v -Vj ) , (17)
which is a diagonal quadratic program with linear constraints, and can be effectively solved via the well-known simplex method, for example.

C Analytic upper bound for regularization parameter in sparse-coding

Let V be the current dictionary at time t ≥ 0 (with the t subscript dropped for ease of notation), and consider the sparse-coding problem

u t ← argmin u∈R k 1 2 Vu -X t 2 2 + α u 1 , (18)
which is equivalent to problem (4), with the choice of penalty φ = . 1 on the codes. Now, it follows from the Lasso theory that u t above will be the zero vector if α ≥ α max (t), where

α max (t) := V T X t ∞ = max 1≤j≤p | V j , X t |
Now, it is clear that

α max (t) ≤ max 1≤j≤p V j 1 X t ∞ ≤ X t ∞ ≤ sup t≥0 X t ∞ ,
the first inequality being a consequence of the Cauchy-Schwarz and the second is due to the constraints on the dictionary atoms V

In particular, this means that a cross-validation procedure for selecting α only need to consider values in the range 0 ≤ α ≤ X ∞,∞ , which depends only on the input data matrix X.

  of the overall algorithm. The Convergence of our algorithm (to a local optimum) is guaranteed since all hypotheses of [18] are satisfied. For example, assumption (A) is satisfied because fMRI data are naturally compactly supported. Assumption (C) is satisfied since the ridge-regression problem (4) has a unique solution. More details are provided in Appendix A.

Figure 1 :

 1 Figure 1: Influence of model parameters. In the experiments, α was chosen according to (9). Left: Percentage explained variance of the decomposition, measured on left-out data split. Right: Average normalized sparsity of the dictionary atoms.

  Mean explained variance of the different models on both training data and test (left-out) data. N.B.: Bold bars represent performance on test set while faint bars in the background represent performance on train set. Pi ct ur e vo ca b. En gl ish re ad in g Pe nn . M at rix Te st St re ng th En du ra nc e Pi ct ur e se q. m em . D ex te rit y Predicting behavioral variables of the HCP[11] dataset using subject-level Z-maps. N.B.: Bold bars represent performance on test set while faint bars in the background represent performance on train set.

Figure 2 :

 2 Figure 2: Main results. Benchmarking our proposed Smooth-SODL (2) model against competing state-of-the-art methods like SODL (sparse online dictionary-learning)[START_REF]Online learning for matrix factorization and sparse coding[END_REF] and CanICA[START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF].

  Vj ) * (I -γA t [j, j] -1 Λ)(v -Vj ) = F * argmin v∈C p , Vj ) * (I -γA t [j, j] -1 Λ)(v -Vj ) = F * P E ( Vj )

j 1 ≤ 1

 11 then u t = 0, ∀t ≥ 0.

  convex subset of R p like an 1 -ball B p, 1 (τ ) or a simplex S p (τ ), defined respectively asB p, 1 (τ ) := {v ∈ R p s.t |v 1 | + |v 2 | + . . . + |v p | ≤ τ } , and S p (τ ) := B p, 1 (τ ) ∩ R p + .Other choices (e.g ElasticNet ball) are of course possible. The radius parameter τ > 0 controls the amount of sparsity: smaller values lead to sparser atoms.• Finally, Ω Lap is the 3D Laplacian regularization functional defined by

Table 1 :

 1 Learning-curve for boost in explained variance of our proposed Smooth-SODL model over the reference SODL model. Note the reduction in the explained variance gain as more data are pooled.

	Smooth-SODL (2) gain factor

For example, see[START_REF] Dohmatob | Benchmarking solvers for TV-l1 least-squares and logistic regression in brain imaging[END_REF][START_REF] Varoquaux | FAASTA: A fast solver for total-variation regularization of ill-conditioned problems with application to brain imaging[END_REF], implemented as part of the Nilearn open-source library Python library[START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF].

Acknowledgment. This work has been funded by EU FP7/2007-2013 under grant agreement no. 604102, Human Brain Project (HBP) and the iConnectome Digiteo. We would also like to thank the Human Connectome Projection for making their wonderful data publicly available.


Cyril Jean

Laurent Belliard

Thomas W Cornelius

Olivier Thomas

Yan Pennec

Marco Cassinelli

Maria Eugenia Toimil-Molares

Bernard Perrin

Spatiotemporal Imaging of the Acoustic Field Emitted by a Single Copper Nanowire

The monochromatic and geometrically anisotropic acoustic field generated by 400 nm and 120 nm diameter copper nanowires simply dropped on a 10 µm silicon membrane is investigated in transmission using three-dimensional time-resolved femtosecond pump-probe experiments. Two pump-probe time-resolved experiments are carried out at the same time on both side of the silicon substrate. In reflection, the first radial breathing mode of the nanowire is excited and detected. In transmission, the longitudinal and shear waves are observed. The longitudinal signal is followed by a monochromatic component associated with the relaxation of the nanowire's first radial breathing mode. Finite Difference Time Domain (FDTD) simulations are performed and accurately reproduce the diffracted field. A shape anisotropy resulting from the large aspect ratio of the nanowire is detected in the acoustic field. The orientation of the underlying nanowires is thus acoustically deduced.

Coherent acoustic phonons with very high in-depth spatial resolution can be generated through the absorption of a femtosecond laser pulse by a nano-structure, such as a semiconducting or metallic thin film [1] or multiple quantum wells [2]. Such generated acoustic waves exhibit a wavelength in the order of magnitude of tens of nanometers, which is proved useful for non-invasive ultrasonic imaging in biological tissues [3,4] or for very high frequency optomechanics [5]. Nonlinear effects could even increase this in-depth spatial resolution by decreasing the spatial extent of the acoustic pulse [6]. However, the lateral resolution of such laser ultrasonic measurements is still limited by the diffraction of light that imposes a micrometric lateral resolution. In order to improve the lateral resolution down to the sub-micrometer scale, various approaches have been suggested. For example, near-field optical miscroscopy (NSOM) was incorporated into the picosecond ultrasonics experimental setup [7,8]. In addition, an opto-acoustic scanning acoustic microscope has also been reported, claiming a lateral resolution below 100 nm [START_REF] Maris | Enhanced ultra-high resolution acoustic microscope[END_REF][START_REF] Maris | Opto-acoustic methods and apparatus for performing high resolution acoustic imaging and other sample probing and modification operations[END_REF]. A Stimulated Emission Depletion (STED) inspired approach has also reduced the acoustic spot lateral dimension to 140 nm [11]. The use of nanostructures as acoustic transducers intergrated in scanning acoustic microscopy systems is also promising. Yang et al. [12] reported the use of a nanowire as an acoustic transducer to optically generate and guide an acoustic pulse. This nanowire can be integrated on an AFM cantilever used to sweep the underlying sample. However, although confined acoustic modes in individual metallic nano-objects is a well-studied topic [13][14][15][16][17], the propagation of acoustic waves in nanometric waveguides [18][19][20][21] and the charac- * cyril.jean@insp.upmc.fr terization of the ability to generate an acoustic wave into the substrate is scarce [22].

In this letter, we propose and demonstrate the use of copper nanowires as efficient nanometric, monochromatic, tunable and anisotropic acoustic transducers. First, investigations of the radial breathing mode of the nanowires of 120 nm and 400 nm diameter are performed using time-resolved pump-probe spectroscopy in the common front side reflectivity setup. Then, a back side transmission setup is used to detect any acoustic feature generated by the nanowire and reaching the other side of the silicon membrane. This approach based on the spatial imaging of the acoustic field enables us to highlight the shape anisotropy of the nanowire-generated acoustic field. Finite Difference Time Domain (FDTD) simulations are compared to the experimental acoustic field data and may allow a rough estimation of the nanowiresubstrate contact size. The imaging of the transmitted 8 GHz and 27 GHz monochromatic acoustic waves through a 10 µm silicon membrane is a very important step towards the use of nano-objects as acoustic transducers for acoustic imaging with nanometric lateral resolution.

Because of their high aspect ratio, their quasi-perfect cylindricity and their smooth contour surface, we have chosen copper nanowires grown by electrodeposition in polycarbonate etched ion-track membranes [23] as model systems to study the confined [24] and propagating [18] acoustic phonons in single nano-objects. The details of their fabrication are described in the supplementary information. Here, these copper nanowires are considered as monochromatic acoustic transducers when simply in contact with a silicon substrate. Indeed, as previously shown [24][25][26], the huge increase in the quality factor of a nanowire's radial breathing mode when suspended in air is interpreted as the suppression of the dissipation channel through the nanowire-substrate contact. In this work, we go one step further, and are able to detect the dissipated acoustic energy in transmission after its propagation through the substrate. We investigate two types of nanowires with length L > 10 µm and diameter d = 120 and d = 400 nm (Figure 1a). According to the continuum mechanics model that holds at these sizes [27,28], the expected frequencies f of our mechanical nanoresonators are in the range f ∼ c L /d ∼ 10 GHz to 30 GHz. At room temperature, these high frequency acoustic phonons suffer severe attenuation in silicon. Consequently the experiments are performed on 10 µm thick silicon membranes obtained by a wet etching in KOH on a 2 inches silicon on insulator (SOI) wafer (Figure 1b,c andd). Experimental details on the fabrication are provided in the supplementary information (Figure S3).

The elastic response of the 400 nm diameter copper nanowires is investigated using a femtosecond optical pump-probe technique described extensively in previous works [29] and in the supplementary information. The mechanical vibrations of the copper nanowires are first studied using the common front side setup where both the pump and the probe laser beams are colinear and focused using the same objective (Figure 2a andS4). Using preliminary SEM experiments (Figure 1b), a single nanowire is located and the subsequent single particle experiment is performed. In this approach, a first blue pump pulse excites the conduction electrons of the metallic nano-object. The subsequent electron cooling by electron-lattice interactions [30] results in the heating of the crystal lattice finally launching the dilatation of the nanowire [31,32]. Because of the electrons fast diffusion, this ultrafast heating is not confined to the surface of the nano-object, the whole width of the nanowire subsequently undergoes periodic radial contrac- tions and expansions. These isotropic volume oscillations correspond to the fundamental breathing mode of the nanowire [24,26]. These periodic lattice oscillations induce a modulation of the metal refractive index, which shows up by a periodic variation of the relative probe reflectivity ∆R/R (Figure 2b). The oscillatory signal is superimposed to an exponentially decaying background as the nanowire cools down by transferring its energy to the substrate. The oscillation frequency f 0 = 8.0 ± 0.2 GHz is extracted from the time trace by either a fast Fourier transform (FFT) or a numerical fitting with a damped sinusoidal function superimposed to an exponentially decaying background. Both methods lead to a quality factor Q on the order of 20. This frequency is compared to the theoretical first radial breathing mode frequency of a free, elastically isotropic and infinite copper cylinder,

f 0 = φ 0 v L /(2πa)
, where a = 200 ± 15 nm is the radius of the nanowire and φ 0 is the first root of the equation ( 1) x (µm)

Time (ns) f) e) d) c) b) a) r(t) 1 t L1 t S t L1 t S t L1 t = 1.32 ns
v T v L 2 J 1 (φ n ) = φ n J 0 (φ n ) (1)
where

v L = 4.5 × 10 3 m • s -1 and v T = 2.2 × 10 3 m • s -1
are respectively the longitudinal and transversal speed of sound in bulk polycrystalline copper [33]. J 0 and J 1 are respectively the Bessel functions of the first kind of order 0 and 1. The theoretical breathing mode frequency reaches f 0 = 8.0 ± 0.3 GHz, where the frequency uncertainty is dictated by the 15 nm fluctuations in the nanowire diameter. Owing to the radial displacement of this breathing mode, the acoustic energy of the nanowire is expected to leak towards the substrate mainly as a longitudinal wave. By flipping a mirror, the experimental setup is thus switched to the backside configuration (Figure 2c andS4). The nanowire remains excited by the blue pump laser, whereas the red probe laser is now focused on the other side of the silicon membrane. The relative reflectivity trace recorded using a Michelson interferometer is plotted in Figure 2d. The origin of the plot corresponds to the electronic coincidence in the frontside configuration. The longitudinal deformation at t L1 = 1.18±0.02 ns is followed by a monochromatic oscillation at the same frequency f 0 = 8.0 ± 0.4 GHz as observed in the frontside configuration (inset Figure 2b). As expected, these oscillations present a quality factor close to that observed in the frontside configuration even if an accurate comparison is not reachable because the damping time of the oscillations is on the order of magnitude of the time delay between the two longitudinal echoes. The second longitudinal echo appears at t L2 = 3.56 ± 0.02 ns after an additional round trip. This second echo is also followed by oscillations. Owing to the well known longitudinal speed of sound along the < 100 > direction of silicon, v SI L = 8.43 × 10 3 m • s -1 , one gets the thickness h = (t L2 -t L1 ) v SI L /2 = 10.0 ± 0.2 µm of the silicon membrane. The shear wave is thus expected at a time delay t S = h/v SI T = 1.84±0.02 ns, where v SI T = 5.84×10 3 m • s -1 is the transverse speed of sound along the < 100 > direction of silicon. As the shear wave reaches the surface while the 8.0 GHz longitudinal oscillations are still occuring, the exact time of arrival cannot be accurately determined. But the main reflectivity change produced by the shear wave appears between 1.8 ns and 2.8 ns. This wide shear wave signature will be discussed later on. The perfect agreement between the resonance frequency of the nanowires' radial breathing mode and the frequency of the generated monochromatic longitudinal wave across the membrane is strong evidence that the nanowire acts as a nanoacoustic wave transducer. We noticed this behavior on several different nanowires : whenever the front side configuration signal shows an oscillatory component, the reported emission phenomenon is observed.

A control experiment is additionally performed. Starting with the time-resolved interferometric measurement in Figure 2d (reproduced in the Figure 3f) when the blue pump beam hits the copper nanowire, the sample is smoothly moved laterally so as the pump beam does not excite the nanowire anymore (Figure 3c). In this configuration, the monochromatic oscillation disappears and the three remaining dipolar peaks at t L1 , t L2 and t S correspond to the longitudinal (direct and after one round trip) and transverse echoes. Both longitudinal and shear waves are effectively excited in the (100) oriented silicon membrane due to the small size of the pump and probe beams [34,35]. The amplitude of the acoustic echoes is similar between the two experiments but the pump power is 5 mW in Figure 3c and 0.5 mW in Figure 3f. The metallic nanowire is thus accountable for the monochromatic oscillations observed in Figure 3f. The surface displacement produced by the acoustic wave was recorded in time and space by displacing the probe laser which is mounted on a piezoelectric stage allowing a movement of 150 µm with an accuracy of 5 nm. Complementary to the experiments described previously (see Figures 3c andf), two different arrangements of the backside configuration were realized where the pump laser beam excites either the nanowire (Figures 3d,e) or the silicon substrate directly (Figures 3a,b). The first experiment is a two dimensional xy surface displacement mapping at a fixed time delay t map = 1.32 ns after the electronic coincidence in the frontside configuration. Considering the micrometric penetration depth of the light in the silicon, the generation of acoustic waves directly in the silicon substrate is inefficient, explaining the extremely poor contrast in the 2D-spatial mapping shown in Figure 3a (the dotted circle is a guide for the eye). On the contrary, when the acoustic wave is generated by a laser-excited copper nanowire, the generation of concentric acoustic waves is highlighted in Figure 3d. These experiments can be performed at different time delays in order to obtain a movie (Movie S1 and S2). As the patterns possess a quasi circular symmetry, another experiment is performed to obtain informations on how the patterns change as the acoustic waves evolve in time. This second experiment is a twodimensional x-t surface displacement mapping at a fixed position y 0 = 7 µm. The alternation of white and black stripes in Figure 3e which is missing on the reference Figure 3c is characteristic of the longitudinal monochromatic acoustic wave emitted by the copper nanowire. In addition, using the rough isotropic approximation, the longitudinal branches in the Figure3 b and e should follow the equation r(t) = ((v L (tt L1 ) + h) 2 -h 2 ) (dotted line in the Figure 3b). This approximation underestimates the experimental r(t). To investigate the nanowire-mediated generation of acoustic waves beyond the rough isotropic estimation, anisotropic FDTD experiments are calculated using a Fortran code developped at the IEMN in Lille. In order to design the geometry of the acoustic source in this simulation, an estimation of the spatial dimensions of the acoustic source generated by the nanowire deposited on the silicon frontside is needed. For now, the nanowires' excitation is produced by a blue light pulse (λ = 400 nm) using a 100× objective with a numerical aperture NA = 0.95. The beam diameter at 1/e 2 is thus measured to be slightly smaller than 1 µm. The excited surface of the silicon substrate may be estimated as the surface of contact between the excited area of the 400 nm nanowire and the silicon. Here, the acoustic source along the nanowire axis is confined to the irradiated part. Possible enlargement associated to guided acoustic phonons [18] along the axis could be neglected due to the huge elastic coupling with the substrate. In a first rough approximation, this source is described in the FDTD simulation by imposing a boundary speed equation

v zz (x, y, z, t) = sin(2πf 0 t) exp 8 (y -y 0 ) 2 /σ 2 y + (z -z 0 ) 2 /σ 2 z + x 2 /σ 2 x
where σ x = 1 µm is related to the pump beam diameter, σ y = 0.2 µm is given by the nanowire diameter and σ z = 0.1 µm is chosen to be small but nonzero to avoid abrupt discontinuity. The out of plane surface displacement u z is thus recorded and plotted in Figures 3g-i using the exact same representation as in Figures 3d-f. The time of flight of the longitudinal and transverse echoes (Figure 3i), as well as the concentric rings (Figure 3g) or alternating white and black stripes (Figure 3h) are in good agreement with the experimental data in figures 3d-f. Modelling the copper nanowire as a submicronic monochromatic longitudinal acoustic transducer seems to effectively describe our experimental observations, even though the transverse acoustic wave signature cannot be grasped perfectly by these simulations. According to the simultaneous experiments in the frontside and the backside configurations of Figure 2, the generated acoustic wave is due to the laser-excited radial breathing mode of the copper nanowire. As the first radial breathing mode frequency, f 0 is inversely proportional to the radius a of the nanowire, it is expected that using smaller nanowires will result in a higher frequency acoustic transducer. Figure 4 shows the results obtained using copper nanowires with 120 ± 10 nm diameter. For these thinner nanowires, the radial breathing mode frequency amounts f 0 = 26.5 ± 0.2 GHz (Figure 4d), which corresponds to a temporal period T 0 = 38 ± 3 ps. Similar to Figure 2, both reflectivity frontside (Figure 4a) and interferometric backside measurements (Figure 4b) are shown. Once again, the monochromatic oscillations obtained using the frontside optical setup (Figure 4 relative reflectivity time trace a) or the backside optical setup (Figure 4 relative reflectivity time trace b) appear to be at the same frequency f 0 = 26.5 ± 0.2 GHz (Figure 4d). Another oscillatory component at ∼ 3 GHz starts around 1.8 ns, at the time delay where the shear wave is expected. This low frequency modulation was also present in Figure 2 at the same time delay but at a lower frequency around ∼ 1.5 GHz. This could be the signature of a low frequency shear acoustic mode generated by the nanowire as the frequency depends on the nanowire's di-ameter. However, as this acoustic feature is missing in the frontside experiment, we are not able to attribute this behavior to any specific acoustic mode. It could be related to some contact mode in the same manner that was previously observed on single gold nanoparticles [36]. Thus, further investigations are still needed to address this question. The insets in these two sub-figures are the derived signals to get rid of the low frequency features in the relative reflectivity signals. Unlike the 8 GHz case, these insets permit a direct and convincing comparison of the damping time of the oscillations between the two cases as the oscillation lifetime is shorter than the time delay between the two longitudinal echoes. Thus, this damping time can be accurately determined and reaches τ = 0.2 ns in transmission to be compared to τ = 0.3 ns in the frontside configuration, which corresponds to a quality factor Q = πf 0 τ ∼ 30 in good adequation with previous studies on supported metallic nanowires [24,37]. A pseudocolor xt spatiotemporal surface displacement mapping along the spatial direction x is presented in the Figure 4c. The white and black stripe alternance, which corresponds to the nanowire breathing mode leaking through the substrate, is only distinguishible directly in front of the generated acoustic wave. As the probe is moved along the x axis, this acoustic feature disappears under the noise. The square root shape of the branches is only visible for the first (and more intense) oscillation. As the wavelength of the 26.5 GHz acoustic wave reaches ∼ 0.3 µm, a spatial imaging (here along the direction x) becomes highly challenging as the probe is focused on a ∼ 1 µm diameter spot. The attenuation that increases with the square of the frequency is also a critical parameter preventing us to obtain nice spatial mapping images.

In the last part of this article, we focus on the 400 nm diameter nanowires vibrating at 8 GHz. As previously described, the blue pump pulse is focused on a nanowire using a 100× objective with a numerical aperture NA = 0.95. The beam diameter at 1/e 2 is thus estimated to be under 1 µm. Consequently, the contact surface between an excited 400 nm diameter copper nanowire and the silicon substrate is estimated as a 200 nm × 1 µm rectangle. This elongated acoustic source should produce a geometrically anisotropic emission, due to a smaller size compared to the acoustic wavelength emitted in the silicon. Thus, the orientation of the underlying nanowire should be detected using such an anisotropic acoustic field. However, the x-y surface displacement mapping of the acoustic field produced by a 400 nm diameter nanowire in Figure 3d does not seem anisotropic : the inner ring seems circular and the color along its perimeter is constant and modulated by a strong noise. The FDTD simulation in Figure 3g does not reveal any obvious anisotropy in the acoustic field either. To understand this observation, let us remind the conditions to observe a diffraction phenomenon. Generally speaking, this phenomenom occur if the size of the obstacle or the slit encountered by the wave is close to the spatial wavelength of the wave. In the present case, the upper limit dimensions of the acous- tic source is a 200 nm × 1 µm rectangle. In comparison, the wavelength of a longitudinal acoustic wave at 8 GHz is λ 0 = v L /f 0 ≃ 1 µm. Both dimensions of the acoustic source are thus below the caracteristic length of the acoustic wave. A significant diffraction is therefore expected in both directions. The remote position of the observer is a second necessary condition to observe a far field diffracted patern. In our experiments, the diffraction will be significantly stronger if the two dimensional xy surface displacement mapping is undergone at a larger scale than the 15 µm × 15 µm of Figure 3. To do so, a 27 µm × 27 µm xy surface displacement mapping is observed at a time delay t d chosen to be 700 ps after the first longitudinal echo detection in the center of the images (Figures 5a, c ande). In order to avoid any confusion with an elastic anisotropy resulting from the silicon substrate, three different nanowires are investigated. A blue line corresponding to the effective orientation of the long axis of the nanowire is added in the graphs. The orientation of the nanowire is checked using a CCD camera positioned on the probe beam path of the frontside setup. The effective position of the nanowires according to the recorded images is thus confirmed without any ambiguity. If the first one (Figure 5a) is arbitrarily referred to as horizontal (0 • ), the second one is at 45 • and the third one is at -55 • . The three corresponding pseudocolor plots of the 27 µm × 27 µm surface displacement mapping are prensented in the Figure 5a,c and e. A substantial modification is also made on the optical setup. Instead of using a 100× objective with a numerical aperture NA = 0.95 to focus the blue pump laser beam, a 50× objective with a numerical aperture NA = 0.5 is preferred. The beam diameter at 1/e 2 is now estimated to be smaller than 2 µm. The upper bound of the contact surface between the excited part of the nanowire and the silicon membrane is therefore a 200 nm × 2 µm rectangle. Only one dimension of this rectangular source is now smaller than the wavelength λ 0 of the considered acoustic wave. This dimension should be favorable to the observation of the diffraction. On the contrary, regarding the 2 µm direction of the source, the phenomenon should be missing. While in Figure 3d when the laser beam diameter is slightly under 1 µm we observe concentric rings, the images in Figures 5a, c and e exhibit incomplete rings. Being confident in the nanowires' orientation, the contrast between each ring appears to be maximal in the direction symbolized by an orange line and minimal in the orthogonal direction symbolized by the blue line corresponding to the long axis of the nanowire. These observations are in good agreement with the fact that diffraction is far more efficient in the direction perpendicular to the nanowire. Moreover, the two profile plots related to each image are depicted on the right hand side of the Figure 5. The spatial period between each oscillation decreases as the distance from the epicenter increases. Using the same isotropic approximation as above, it can be shown that the the spatial separation ∆R between two successive rings is nonlinear (see the supplementary information for details). In the case investigated here, the function decreases from 2.4 to 1.5 µm from the inner to the outer rings. However, the radius of the rings are underestimated (3.9 µm, 6.2 µm, 8.0 µm and 9.6 µm in comparison to the experimental values 6.0 ± 0.2 µm, 8.7 ± 0.3 µm, 10.1 ± 0.3 µm and 11.7 ± 0.3 µm). If this isotropic estimation grasps the good tendency, an anisotropic model is required for quantitative analysis. The same FDTD simulations are performed, modifying only σ x to σ x = 2 µm compared to Figure 3g in the above study. The simulated pseudocolor plot and the corresponding profiles are displayed in Figure 5g and h, respectively. The calculated radiuses of the different rings amount to 6.0 µm, 8.5 µm, 10.3 µm and 11.7 µm, in excellent agreement with the experimental radiuses. The spatial separation effectively decreases between two successive rings from 2.5 µm, to 1.8 µm and finally to 1.4 µm. The FDTD simulations clearly evidence the elastic anisotropy by the fact the concentric rings structures are not perfectly circular. It has to be differentiated from the shape anisotropy resulting from the elongated source. As a last physical insight, the shape anisotropy is quantified in silico. The simulation is undergone with σ y varying from 100 nm to 2 µm. A plot similar to Figure 5h is obtained and the ratio between the oscillation amplitude along the orange line and the blue line is extracted. As this ratio can be in principle extracted from experimental images, it could be a possible way to estimate σ y and thus to give some insight on the surface contact between the nanowire and substrate which is a challenging field of study [38]. However, below 1 µm, a variation of the surface of adhesion of ∼ 100 nm leads to a relative variation of the ratio inferior to 1% which is not distinguishable owing to the level of noise in our experiments.

To conclude, 400 nm and 120 nm diameter electrochemically grown copper nanowires have been used to generate 8 to 27 GHz, monochromatic and longitudinal acoustic phonons with sub-micronic lateral dimensions. Following a standard time-resolved pump-probe spectroscopy experiment in reflectivity to assign the acoustic frequencies to the radial breathing mode of the nanowires, picosecond laser ultrasonic experiments have been performed in transmission on silicon membranes. The as-generated acoustic waves produced by the vibrating nanowires have been imaged in transmission in both time and space domains. The transmitted acoustic field exhibits a strong geometric anisotropy due to the high aspect ratio of the nanowires. The propagation of these gigahertz acoustic phonons generated by the copper nanowires has been modeled using FDTD simulations which gives a good agreement with the experimental data. The transmission of the acoustic waves of the nanowire through a 10 µm silicon substrate definitely shows the good adhesion of nanometric objects simply dropped on their substrate, producing a strong dissipative channel for the nanowire mechanical oscillator. Furthermore, this technique may be used to quantify the surface contact between nanowires and silicon substrates if stronger signal to noise ratios are achieved. Further investigations could also imply the modification of the substrate-wire coupling in order to exalt or inhibit the reported tunable emission of acoustic waves by a nanowire. This demonstration of nanometric longitudinal wave acoustic transducers could also be useful for improved lateral resolution of acoustic microscopy using for instance synthetic aperture focusing techniques [39]. The implementation of such a nanometric transducer on a sweeping device is now needed to pave the way to acous-tic microscopy with nanometric lateral resolution.

FIG. 1 .

 1 FIG. 1. (a) SEM top view of the etched silicon on insulator wafer. (b) Close SEM view of a 120 nm copper nanowire supported by the silicon substrate. (c) 3D schematic of the etched 100 µm SOI wafer, resulting in a 10 µm silicon membrane with a squared surface topped with a thin 100 nm aluminium layer. (d) Side schematic view of a nanowire resting at the bottom of a silicon membrane.

FIG. 2 .

 2 FIG. 2. (a) Front side configuration. Both the blue pump laser and the red probe laser are focused on a single 400 nm diameter copper nanowire using the same 100× objective with a numerical aperture NA = 0.95. (b) Relative reflectivity ∆R/R obtained in the front side configuration Figure 1a. It exhibits an oscillatory signal with a period T = 125 ps. Inset : normalized FFT power spectrum of the derived signal of Figure 2d (in red) and Figure2b (in green). (c) Back side configuration. The two laser beams are focused on different sides of the silicon membrane. The probe is mounted on a piezoelectric stage in order to map the displacement field on a maximum 150 µm × 150 µm area. (d) Relative reflectivity ∆R/R obtained in the back side configuration Figure 1c. It exhibits an oscillatory component with a period T = 125 ps.

FIG. 3 .

 3 FIG. 3. (a), (b), (c) Reference mapping : the pump hits directly the silicon membrane (schematic inset Figure 3c). (d), (e), (f) The pump excites the first radial breathing mode of the 400 nm diameter copper nanowire (schematic inset Figure 3f) which generates a monochromatic longitudinal acoustic wave leading to concentric rings. (g), (h), (i) FDTD simulation : the acoustic source is a 8 GHz monochromatic longitudinal displacement imposed on the surface of the silicon layer. (j) FFT power spectrum of the red part of the Figure 3f and the turquoise part of the Figure 3i. (a), (d), (g) Two-dimensional xy surface displacement mapping at tmap = 1.32 ns. (b), (e), (h) Two-dimensional xt surface displacement mapping at fixed y0 = 7 µm. (c), (f), (i) One-dimensional (time domain) surface displacement recording on the position x0 = 10 µm and y0 = 7 µm. The red (or turquoise) part of the signal matches the corresponding red (or turquoise) line on the Figure 3b, e and h.

2

 2

FIG. 4 .

 4 FIG. 4. (a) Relative reflectivity ∆R/R obtained in the front side configuration (See schematic in the inset). Both the blue pump and the red probe are focused on a single 120 nm diameter copper nanowire using the same 100× objective with a numerical aperture NA = 0.95. The other inset is the derived signal exhibiting a 35 ps oscillatory signal. (b) Relative reflectivity ∆R/R obtained in the back side configuration (see schematic in the inset) when the blue pump is focused on a single 120 nm diameter copper nanowire. The other inset is the derived signal exhibiting an oscillatory signal with a periodicity T = 38 ps (c) Twodimensional (spatiotemporal) surface displacement mapping along the spatial direction x. The red line corresponds to the red part of the oscillatory signal in Figure 4b. (d) Normalised FFT power spectrum of the inset of Figure 4b (in red) and Figure 4a (in green).

1 FIG. 5 .

 15 FIG. 5. (a),(c) and (e) Two-dimensional x-y surface displacement mapping 700 ps after the longitudinal echo detection on the central position with a ×50 objective. The three different underlying nanowires are oriented along the blue cyan line. (g) FDTD two-dimensional xy surface displacement (uz) mapping 700 ps after the longitudinal echo detection on the central position. The acoustic source is ∼ 2 µm wide along the blue cyan line and ∼ 0.2 µm wide along the orange line. (b), (d), (f) and (h) Surface displacement profile, respectively along the orange and the cyan line. The offset is for readability purpose.

ACKNOWLEDGMENTS

MET and MC thank the Deutsche Forschungsgemeinschaft (DFG) for financial support within the priority program SPP 1386.