
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38451944
202314697
10.1073/pnas.2314697121
research-articleResearch Articleapp-mathApplied Mathematics404
Physical Sciences
Applied Mathematics
Correlation-informed ordered dictionary learning for imaging in complex media
Moscoso Miguel a https://orcid.org/0000-0001-8461-1578

Novikov Alexei b
Papanicolaou George papanicolaou@stanford.edu
c 1
Tsogka Chrysoula ctsogka@ucmerced.edu
d 1 https://orcid.org/0000-0002-5839-8889

aDepartment of Mathematics, Universidad Carlos III de Madrid, Leganes, Madrid 28911, Spain
bDepartment of Mathematics, Pennsylvania State University, University Park, PA 16802
cDepartment of Mathematics, Stanford University, Stanford, CA 94305
dDepartment of Applied Mathematics, University of California, Merced, CA 95343
1To whom correspondence may be addressed. Email: papanicolaou@stanford.edu or ctsogka@ucmerced.edu.
Edited by David Weitz, Harvard University, Cambridge, MA; received August 24, 2023; accepted February 4, 2024

7 3 2024
12 3 2024
7 9 2024
121 11 e231469712124 8 2023
04 2 2024
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

We show how to exploit large and diverse sets of data to generate high-quality images in complex media. Previous imaging methods lose resolution in complex media because they are based mostly on homogenous reference media. Our approach produces images whose resolution is comparable to that of a homogeneous medium. We use algorithms that come from computer science and statistics rather than from imaging.

We propose a method for imaging in scattering media when large and diverse datasets are available. It has two steps. Using a dictionary learning algorithm the first step estimates the true Green’s function vectors as columns in an unordered sensing matrix. The array data comes from many sparse sets of sources whose location and strength are not known to us. In the second step, the columns of the estimated sensing matrix are ordered for imaging using the multidimensional scaling algorithm with connectivity information derived from cross-correlations of its columns, as in time reversal. For these two steps to work together, we need data from large arrays of receivers so the columns of the sensing matrix are incoherent for the first step, as well as from sub-arrays so that they are coherent enough to obtain connectivity needed in the second step. Through simulation experiments, we show that the proposed method is able to provide images in complex media whose resolution is that of a homogeneous medium.

imaging
complex media
dictionary learning
multidimensional scaling
Spanish AIE PID2020-115088RB-I00 MIGUEL ANGEL MOSCOSOAlexei NovikovGeorge PapanicolaouChrysoula Tsogka AFOSR FA9550-23-1-0352 MIGUEL ANGEL MOSCOSOAlexei NovikovGeorge PapanicolaouChrysoula Tsogka AFOSR FA9550-23-1-0352 and 443 FA9550-23-1-0523 MIGUEL ANGEL MOSCOSOAlexei NovikovGeorge PapanicolaouChrysoula Tsogka AFOSR FA9550-23-1-0352 and 446 FA9550-21-1-0196 MIGUEL ANGEL MOSCOSOAlexei NovikovGeorge PapanicolaouChrysoula Tsogka
==== Body
pmcHigh-resolution imaging in complex media faces challenges due to wavefront distortion caused by scattering from inhomogeneities. In this paper, we introduce a method for imaging in a randomly inhomogeneous medium involving two basic algorithms. The first is a sparse dictionary learning algorithm (1) in order to estimate Green’s function vectors between focal or source points in the image window and receiver locations on the array. The second is a multidimensional scaling (MDS) algorithm (2) to convert information about correlations of Green’s function vectors into positions of the focal points in the image window.

To accomplish the first step, we use a sparsity-promoting modification of the method of optimal directions (MOD) (3) to learn an (unordered) dictionary of Green’s function vectors that characterize the propagation of signals from a set of focal points, or sources in the image window, to the array. Here, unordered means that we do not know which focal points are associated with the estimated column vectors of the dictionary. In this step, we assume that an abundance of sensing measurements is available. Specifically, we have access to measurements for multiple signals emanating from many sparse sets of sources, but we do not have prior knowledge of their locations or amplitudes. This dictionary learning algorithm enables us to estimate Green’s function vectors with high accuracy under the condition that these vectors are sufficiently incoherent, which means that their normalized inner product is sufficiently small. Given a configuration of sources or focal points in the image window, this implies that the receiver array must be large enough. We present this dictionary learning step in Section 2.1.

The goal of the second step is to associate each Green’s function vector with its corresponding focal point in the image window, which means that we want to find the correct order of the columns in the estimated matrix of Green’s function vectors. We could back-propagate these vectors into the image window using a reference homogeneous medium. This is the Kirchhoff’s migration approach that only works well when the fluctuations are weak (4). For a given set of Green’s function vectors, we could also try to estimate the position of their focal points using source localization algorithms, commonly used in wireless communications (5–8). However, these algorithms use information based on distances between sources and receivers and are therefore very sensitive to noise. They are not suitable for imaging in media with strong fluctuations.

Instead of estimating the distance between a focal point and a receiver, we can obtain a much more accurate estimate of the distance between two nearby focal points. The correlation of the estimated Green’s function vectors gives such an estimate but, of course, we do not know where their focal points are located in the image window. By cross-correlating each Green’s function vector with all the others, we can identify its nearest neighbors. This is a key observation that allows us to generate a proxy distance between column or Green’s function vectors by counting the smallest number of neighborhoods that connect them. This provides a connectivity-based proxy distance between all pairs of column vectors that we can use with the multidimensional scaling (MDS) algorithm (2) for identifying Green’s function vectors with their focal points up to a rotation, translation, and scaling. The resulting relative configuration of points can be spatially fixed with a few (two or three) known reference points in the two-dimensional image window. The use of a proxy metric based on connectivity is done by the MDS-MAP algorithm (9, 10). Constructing the connectivity-based proxy distance using cross-correlations is described in Section 2.2. We should point out that this use of cross-correlations is different from another use of cross-correlations in imaging, in fields such as geophysics (11, 12), where the gathered data are cross-correlated to estimate the Green’s function between two receivers in complex media.

1. Imaging Problem Setup

Suppose that an array of N receivers records waves generated by sources located over a region of interest, called the image window. The receivers are located at points r→j, and the sources at unknown locations z→i. In Fig. 1, it is assumed that the array is one-dimensional. The coordinates parallel to this array are the cross-range coordinates, and the ones orthogonal to it are the range coordinates. The medium between the array and the unknown sources fluctuates randomly in space as illustrated in Fig. 1. The Green’s function that characterizes wave propagation of a signal of frequency ω from a point z→ to a point r→ in the random medium satisfies the wave equation[1] ΔG(z→,r→)+κ2n2(z→)G(z→,r→)=δ(z→−r→),

Fig. 1. Schematic of the data collected on the array when three sources located at z→i, z→i′ and z→i″ simultaneously emit signals. It is the superposition of these signals that are recorded on the array of receivers located at points r→j, j=1,…,N.

where κ=ω/c0 is the wavenumber with constant reference wave speed c0. The random index of refraction is n(z→)=c0/c(z→) with local wave speed c(z→). In a homogeneous medium, c(z→)≡c0 for any location z→ and, in this case, G(z→,r→)=G0(z→,r→), where[2] G0(z→,r→)=exp(iκ|z→−r→|)4π|z→−r→|.

In random media, however, the wave speed c(z→) depends on the position z→. We consider a variable wave speed satisfying[3] 1c2(z→)=1c02(1+σμ(z→l)),

where l is the correlation length of the inhomogeneities that is characteristic of their size. In Eq. 3, σ determines the strength of the fluctuations around the constant speed c0, and μ(·) is a stationary random process with zero mean and normalized autocorrelation function R(|z→i−z→i′|)=E(μ(z→i)μ(z→i′)), so R(0)=1.

We write the data received on the array of N receivers in vector form, with Green’s function vector[4] g(z→)=[G(r→1,z→),G(r→2,z→),…,G(r→N,z→)]T,

and we introduce the N×K sensing matrix [5] G=[g(z→1)⋯g(z→K)]

defined on a grid {z→i}i=1,⋯,K spanning the image window. The sensing matrix G maps a distribution of sources in the image window to the (single frequency) data received on the array. The multi-frequency case is described in Section 3. For a given configuration of sources on the grid represented by the vector x∈CK, the data recorded on the array is given by [6] y=Gx,

where y∈CN. Here, x is a vector whose jth component represents the complex amplitude of the source at location z→j in the image window, j=1,⋯,K.

Because the medium is random, the sensing matrix G in Eq. 6 is not known. Our imaging problem is to estimate this matrix from a set of M samples or observations {yi}i=1,⋯,M, with yi=Gxi. The number of observations is large with respect to the dimension of the vectors xi, i.e., M≫K. Note that xi is also unknown but we assume that it is sparse, implying that the samples yi can be represented as a linear combination of a small number of columns of the unknown sensing matrix G. Since we do not know xi, both the locations and the amplitudes of the sources are unknown. We assume that for the few sources that are active for every sample the modulus of their amplitude takes values in a bounded interval away from zero.

We note that, for imaging problems, the coherence between the columns of the sensing matrix G increases as the grid in the image window becomes finer. This can be challenging for the sparse dictionary learning algorithm, described in the next section, as its convergence is guaranteed under incoherence or restricted isometry property assumptions (13). Coherence is defined as the maximum of the normalized inner product between different columns of the matrix, i.e.,[7] ν= maxi≠ji,j=1Kg(z→i)∗g(z→j)||g(z→i)||  ||g(z→j)||.

Given the size of the recording array and the bandwidth, in order to limit the coherence between the columns of G, we assume that the grid in the image window is not finer than the support of the point spread function in a reference homogeneous medium; see Fig. 2A. This is essential for the first step of the imaging method, the dictionary learning step. The point spread function of an imaging system is the image one obtains when the signal from a single source is used as input. The proposed two-step imaging method is described next.

Fig. 2. Refocused spots with (A) a large array used in the first step of the algorithm and (B) a small array used in the second step.

2. The Two-Step Imaging Method

The first step of our imaging method uses a dictionary learning algorithm and allows us to recover the columns of the sensing matrix up to a permutation. Although in most of the applications of dictionary learning the order of the columns is not an issue, it is essential in the imaging case. That is because even though we recover the Green’s function vectors, the imaging problem is still not solved as we do not know their correspondence with the grid points in the image window. To create an image, we need to associate each column of the estimated matrix G^ with the corresponding focal point in the image window. This is challenging since the propagation medium is unknown, so we cannot just back-propagate the recovered Green’s function vectors, as it is done in time reversal.

In the second step of our method, we deal with the focal spot localization problem where we associate each recovered Green’s function vector g^i with its corresponding focal or source point in the image window. The key idea is the use of cross-correlations between the columns of the estimated matrix G^ to identify the nearest neighbors of their focal or source points in the image window and then infer their associated focal points on the grid. Coherence becomes crucial in this step, as we cannot identify the nearest neighbors of each focal point if G^ is column-incoherent. We can do this, however, by using for the cross-correlations a suitable subset of the elements of the estimated columns g^i corresponding to a fraction of the size of the recording array. An example is illustrated in Fig. 2B where a subarray of half the size is used. The grid reconstruction problem is then solved using the MDS algorithm with a proxy distance, as in sensor network localization problems (10). No Euclidian distance information is known but a proxy distance is obtained from connectivity information at limited range that is recovered from cross-correlations.

The cross-correlations formed can be interpreted as time reversal experiments. That is, the signals recorded on the array are re-emitted in the same medium and given the time-reversibility of the wave equation, those signals will focus on the location from which the original signal was emitted. Our connectivity reconstruction relies on this fundamental property of the wave equation and therefore is robust to the complexity of the medium. Our numerical simulations confirm that this is the case.

2.1. Step One: Dictionary Learning.

We discuss here an algorithm aimed at learning from the gathered data a dictionary A∈CN×K that represents the normalized sensing matrix (5). We assume that the recorded signals yi∈CN, i=1,⋯,M, the data, come from only a few sources and can therefore be represented as a linear combination of a small number of columns in the dictionary A we want to determine. This means that yi=Axi, where xi∈CK are sparse vectors that represent unknown collections of sources firing at the same time. For imaging applications, we can assume that the columns of A have unit lengths. We also assume that we know the dimension K (with usually K≫N), which is the number of points in the image window and therefore specifies the resolution of the image. An estimate of K can be based on the resolution of the imaging setup expected in a homogeneous medium. In a random medium resolution in time reversal, but not in imaging, will improve (14, 15) and therefore K could be larger. We assume here that K is chosen based on a homogeneous medium.

To find the dictionary A and the layout of sources, we define the matrix X=[x1,⋯,xM]∈CK×M and the data matrix Y=[y1,⋯,yM]∈CN×M, and solve the problem[8] minA,X‖AX−Y‖F2s.t.‖xi‖0≤s,i=1,⋯,M,

where ‖·‖0 counts the number of nonzero elements and s is the expected sparsity level. The decomposition Y=AX is unique up to permutations of the columns of A and rows of X provided that the data Y is rich enough (16).

How much data do we need, that is, how big should M be? As already noted, in imaging, we usually have K≫N. If the sparsity s is fixed independent of N then the condition[9] M>KlogK

is sufficient for a suitable probabilistic model of A,X,Y (17, 18).

Problem (8) is nonconvex, as the constraint is not convex and both A and X are unknown. However, its solution can be found efficiently by means of an alternating optimization procedure that uses the ℓ1-norm instead of the sparsity count, provided the initialization is close enough to the true solution and the columns of X are sparse enough. We refer the reader to ref. 13 for local convergence guarantees.

Specifically, if A is known, then X can be obtained from the ℓ1-norm minimization problem[10] min‖X‖1  subject to AX=Y,

that can be easily solved by several different algorithms (19–22). Here, we solve Eq. 10 using a Generalized Lagrangian Multiplier Algorithm (GeLMA) (23). In the next step, if X is known, the minimization problem for A[11] minA‖AX−Y‖F2,

can be easily solved. The exact solution is A=YXT(XXT)−1, provided XXT is invertible (actually stably invertible). Once A has been computed, we normalize its columns to one.

To summarize, in order to solve Eq. 8, we alternate between problems Eq. 10 and Eq. 11 to update X and A sequentially one after the other. Each iteration has two steps, starting from an initial guess X0 and A0. At the beginning of iteration l≥1, we have Al−1 and Xl−1, and we solve Eq. 10 with A=Al−1 to obtain Xl using GeLMA. Then, we solve Eq. 11 with X=Xl fixed to find[12] Al=YXlT(XlXlT)−1

in the second step. This is very much like the MOD algorithm proposed in ref. 3 for signal compression that uses Matching Pursuit for finding Xl instead of GeLMA in the first step.

Our numerical experiments indicate that, given a suitable initialization, this dictionary learning algorithm can construct the matrix of Green’s functions for wave propagation in random media. However, the columns of this matrix, the dictionary, are unordered. They cannot be used for imaging because we do not know the points where they focus in the image window. The next section addresses this problem.

2.2. Step Two: Grid Reconstruction.

We describe now an algorithm for finding the focal spots in the image window from the estimated Green’s function vectors {g^i}i=1K. It is the range-free or connectivity-based sensor localization algorithm (9), analyzed in ref. 10. Our main contribution here is to determine connectivity from cross-correlations of the Green’s function vectors {g^i}i=1K, which now must retain some coherence. However, the dictionary learning algorithm of the previous section requires incoherence, which means that the ν in Eq. 7 for the estimated {g^i}i=1K is small. By using a subset of the components of the {g^i}i=1K, corresponding to a subarray as illustrated in Fig. 2, we increase their coherence. In this section, we use the same notation {g^i}i=1K for their subsampled components.

This increased coherence is essential for the connectivity-based localization algorithm because it allows us to introduce a graph G=(V,E) where the vertex set V={1,2,⋯,K} is associated with the estimated Green’s function vectors {g^i}i=1K. A pair of vertices (i,j) is then connected by an edge in E and assigned the value one in the adjacency matrix of the graph, if the cross-correlation g^i∗g^j is sufficiently close to one in absolute value. Otherwise, the pair (i,j) is not connected, and zero is entered in the adjacency matrix. The size of the subarray of receivers is adjusted so that each vertex has up to k=2r edges, where r=2,3 is the ambient dimension of the image window.

The proxy distance between two Green’s function vectors is now the geodesic graph distance between their corresponding vertices. That is, the proxy distance between g^i and g^j, denoted by d^ij, is the number of edges in the shortest path connecting i and j. We use this proxy distance as a replacement of the Euclidean distance between pairs of focal points in the image window associated with Green’s function vectors in the MDS algorithm. The resulting configuration of focal points Z^=[z→^1,z→^2,⋯,z→^K]T in the image window provides an estimate for the true configuration of focal points Z=[z→1,z→2,⋯,z→K]T, up to rotation, translation and scaling. This is the MDS-MAP algorithm (9) with our correlation-based proxy distance as in Algorithm 1.

When the true Euclidean distance D=(dij) is used instead of D^=(d^il), the classical metric MDS algorithm (2) recovers the configuration of focal points Z=[z→1,z→2,⋯,z→K]T up to rotation and translation. In this case, the input is a K×K squared distance matrix D with entries dij, dij=(z→i−z→j)T(z→i−z→j), and the output the K×r configuration matrix of focal points Z. We have that (2)

[13] −12LDL=LZZTL,

where L=IK−1lK1lKT/K is a centering matrix, with IK the K×K identity matrix, and 1lK the column vector of all ones. This means that the matrices ZZT and −D/2 are equal when the center of mass of the configuration is moved to zero. For the Euclidean distance matrix D, Algorithm 1 determines the Euclidean coordinates of the focal points.

The rank of the matrix P=−12LDL equals the ambient dimension r of the image window when the Euclidian distance matrix D is used in the MDS algorithm. When we use the geodesic distance D^ on the graph, then the rank of P is not equal to r any more. However, the first r singular vectors of P are close to the true coordinates Z (up to centering and rotation) (10). This is illustrated in Fig. 3. The absolute location of the focal points in the image window can be determined using the true location of a few of them, the anchors. These anchors allow us to find the proper rigid transformation and scaling to superimpose the given configuration over them. The anchors can be known a priori or their location can be estimated using coherent interferometric imaging (24). The number of anchors needed is small, typically r+1.

Fig. 3. The singular values of the doubly centered distance matrix P normalized by the maximal one are plotted with red circles when D is used (rank is exactly 2) and with blue stars when D^ is used. There are exactly two top singular values in the second case as well, plotted with blue stars, and the lower eigenvalues drop to zero fast but are not immediately zero as with the red circles.

3. Numerical Experiments

To simulate wave propagation in random media, we use the random travel time model (24, 25) and references therein which provides an analytical approximation for the Green’s function in Eq. 1 in the high-frequency regime in random media with weak fluctuations and large correlation lengths ℓ=100λ compared to the central wavelength λ, given by[14] G(z→,r→)=G0(z→,r→)expiσκ|z→−r→|∫01μz→l+sl(r→−z→)ds.

Comparing Eqs. 14 and 2 for a homogeneous medium, we see that, in this regime, only the phases are perturbed by the random medium while the magnitudes remain unchanged. This model is widely used in adaptive optics as it captures well waveform distortions in heterogeneous media (26). It does not capture, however, the delay spread (coda) due to multiple scattering and does not provide cross-range diversity.

In our numerical experiments, the distance between the array and the image window L=100ℓ is large, so the small distortions produced by each inhomogeneity build up over the propagation distance and are significant at the receivers. The strength of the fluctuations σ is scaled by the dimensionless parameter λ/lL, for which the standard deviation of the random phase fluctuations in the Green’s function is O(1). The strength of the fluctuations σ~=σ/(λ/lL) in the simulations is σ~=0.6 or σ~=0.8.

The challenge is to obtain statistically stable results in these media. Some previously proposed methods, such as coherent interferometry, provide stable images but at the expense of resolution, causing some blurring (24). Our approach produces stable images whose resolution is comparable to that of a homogeneous medium.

We consider the following setup for our numerical simulations. In the first step of dictionary learning for the sensing matrix G, we use the multi-frequency data recorded with a large array aperture a=48ℓ with Nr=145 equally spaced receivers (Fig. 2). In the second step for the grid reconstruction, we use the data corresponding to half the array aperture, so a=24ℓ. The bandwidth [0.5f,f], with f=c0/λ, is discretized with Nf=10 equally spaced frequencies, and is the same for both steps of the algorithm. We organize the multiple frequency data column-wise, soY=[Y(f1)⊺,Y(f2)⊺,⋯,Y(fNf)⊺]⊺.

The multi-frequency sensing matrix is nowgi=[g(z→i,f1)⊺,g(z→i,f2)⊺,⋯,g(z→i,fNf)⊺]⊺,

i=1,⋯,K. Thus, the sensing matrix G=[g1⋯gK] has dimensions N×K with N=Nr·Nf. The sampling of the 20×20 points in the image window is based on the homogeneous medium array resolution, O(λL/a) in cross-range and O(c0/B) in range (15, 27). This means that we implicitly assume that the sources are far part in the sense that the random Green’s function vectors corresponding to this discretization are non-coherent, i.e., the coherence ν Eq. 7 is small.

In Fig. 4, we assume that the sensing matrices corresponding to a random medium G and to the homogeneous medium G0 are known, and we show the cross-correlation matrices of G∗G (Left), G0∗G0 (Center), and G∗G0 (Right). For the homogeneous medium, the Green’s function used is given by Eq. 2. Each row i in these images corresponds to a time reversal experiment where a source located at z→i emits a pulse, and the recorded signals are time reversed and emitted back into the medium. When the waves are re-emitted into the same medium in which the measurements were obtained, as in the Left and Center images of Fig. 4, they retrace the original scattering process and arrive back approximately at the point at which they were emitted, that is, the focal point. However, when the back-propagation is done in a different medium, as in the Right image of this figure, there is no re-focusing. In Fig. 4, the large values (lighter blue color) correspond to re-focusing points. This figure shows that a) time reversal of waves into random and homogeneous media are similar and b) that we cannot use the homogeneous medium to recover this structure if there is scattering.

Fig. 4. Cross-correlations of the sensing matrix. Left: cross-correlations of the sensing matrix in the random medium. Center: cross-correlations of the sensing matrix in the homogeneous medium. Right: cross-correlations between the sensing matrix in the random and homogeneous media. Strength of the fluctuations of the random medium σ~=0.8.

In our numerical experiments, we assume that we have a diverse set of data Y=[y1,y2,…,yM], with yi=Gxi. Both the sensing matrix G∈CN×K and the sparse vectors xi∈CK are unknown. We assume data corresponding to a large number of experiments, so M≫K. Given this set of data, we want to recover the columns of the sensing matrix G=[g1⋯gK], whose rank is approximately 200<K=400 and whose coherence is ν=0.7. The sensing matrix is rank-deficient because the resolution of the image window is high, with pixel sizes λL/a in cross-range and c0/B in range.

The results of the first step of the proposed strategy are depicted in Fig. 5 for large (red lines) and small (blue lines) arrays. We solve problem Eq. 8 as described in Section 2.1 and measure the success of this first step as follows. We form, for every (normalized) recovered column g^i, the cross-correlations with all the columns of the true sensing matrix G and represent in Fig. 5 the maximum value[15] Cmax(i)=maxj|g^iTgj|

for a sparsity level s=4 (Left) and s=8 (Right). We observe values very close to 1 in both cases when the columns of the sensing matrices are for large array apertures (red lines). This means that the true Green’s function vectors are recovered when large arrays are used because they are incoherent. However, when smaller arrays are used (blue lines), the Green’s function vectors are coherent and some of them are not recovered. It is important to recover accurately all, or almost all, Green’s function vectors because, otherwise, we cannot establish their connectivity properly and, therefore, we cannot reconstruct the grid in the image window in the second step.

Fig. 5. Maximum correlation between each estimated column and the true ones in the sensing matrix, as in Eq. 15. In red, the results when the columns of the matrix are more incoherent (a=48ℓ). In blue, the results when the columns of the matrix are more coherent (a=24ℓ). Sparsity s=4 on the Left and s=8 on the Right.

Fig. 6 shows the results for the grid reconstructions accomplished with Algorithm 1 described in Section 2.2 using k=4 neighbors. This algorithm provides the correspondence between the Green’s function vectors found in the first step and their focal points in the image window. From Left to Right, we show the results when (Left) all the pairwise Euclidean distances between the focal points are known in a homogeneous medium, (second from the Left) when only Euclidean distances between the four nearest neighbors are known in a homogeneous medium, (second from the Right) using only connectivity information in a random medium with σ~=0.6, and (Right) using only connectivity information in a random medium with σ~=0.8. In all the cases, the sparsity is s=8.

Fig. 6. From Left to Right: Grid reconstruction from true Euclidean distances using the MDS algorithm when all the pairwise distances are assumed known; from true Euclidean distances when only distances corresponding to the four nearest neighbors are assumed known; using the MDS-MAP algorithm with geodesic graph distances for σ~=0.6; and using the MDS-MAP algorithm with geodesic graph distances for σ~=0.8. Sparsity s=8 in all cases.

Algorithm 1 provides grid positions up to a rigid transformation and scaling. In Fig. 7, we post-process the results shown in the second from the right and right images in Fig. 6 to transform them to absolute positions using three anchors. We observe that the grids are quite well reconstructed near the center but bent toward the edges. This occurs because our geodesic graph distance is the scaled l1 distance on the grid, and an embedding of such distances into Euclidean spaces leads to such distortions. Naturally, there is no grid deformation shown in the left image of Fig. 6 since Euclidean distances are used for all focal points.

Fig. 7. Using three points as anchors, i.e., assuming the location of those three points is known, we can estimate the scaling and the rotation needed to recover the absolute grid positions. We compare the recovered locations (red stars) with the true ones (blue circles) where σ~=0.6 (Left) and σ~=0.8 (Right). Sparsity s=8 in both cases.

After the two steps of the proposed strategy, we recover the ordered sensing matrix G^, so we can image any signal measured at the array into the image window. In Fig. 8, we back-propagate a signal g(z→j) from a source located at z→j using the recovered Green’s function vectors. Thus, the image formed at points z→i,i=1,…,K, is[16] I(z→i;z→j)=g^(z→i)∗g(z→j).

Fig. 8. From Left to Right, image formed with Eq. 16 using the true random Green’s functions, the homogeneous Green’s functions and the recovered ones with the proposed method. Here, the sparsity is s=8. The strength of the fluctuations is σ~=0.6 for the Top row and σ~=0.8 for the Bottom row.

As before, the hat in Eq. 16 denotes the recovered Green’s function vectors of the sensing matrix using the two-step method introduced here. As illustrated in Fig. 8, the produced image (Right) is similar to the one obtained using the true Green’s functions (Left) and significantly better than the one obtained using the homogeneous Green’s function (Center). This figure shows the need for recovering accurate estimates of the Green’s function vectors for imaging in random media since the ones corresponding to a reference homogeneous medium provide very noisy, useless images with Eq. 16. The top and bottom rows are images obtained in different realizations of the random media with σ~=0.6 and σ~=0.8, respectively.

4. Discussion

We propose here a data-driven approach for imaging in random media. The data are multiple signals recorded by an array of receivers and coming from many sets of sparse sources whose location and strengths are unknown. Imaging is done with a two-step method. In the first step, we show that dictionary learning algorithms can estimate the random Green’s function vectors as columns of an unordered matrix. Ordering is not an issue in other applications of dictionary learning but it is essential for imaging as we need to know the correspondence between them and the grid points in the image window. This is done in the second step using a multidimensional scaling algorithm that uses local information at each grid point to compute the global layout of the grid points. A key element of the proposed approach is that the local information is obtained from the cross-correlations of the estimated Green’s function vectors. These cross-correlations determine the connectivity between the grid points.

The numerical experiments show that using these two basic algorithms, dictionary learning and multidimensional scaling, we can form stable images in random media without loss of resolution. This demonstration of principle can be extended to imaging in more complex ambient media, but this may require additional steps in the imaging method.

5. Materials and Methods

All the imaging data used in the numerical experiments are generated numerically as described in detail in Section 3. The implementation of the two-step imaging method that we are proposing is also described in detail in this section.

Miguel Moscoso’s work was supported by the Spanish AEI grant PID2020-115088RB-I00. Alexei Novikov’s work was partially supported by AFOSR FA9550-23-1-0352 and FA9550-23-1-0523. The work of George Papanicolaou was partially supported by AFOSR FA9550-23-1-0352. The work of Chrysoula Tsogka was partially supported by AFOSR FA9550-23-1-0352 and FA9550-21-1-0196.

Author contributions

M.M., A.N., G.P., and C.T. performed research and wrote the paper.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

There are no data underlying this work.

This article is a PNAS Direct Submission.
==== Refs
1 K. Kreutz-Delgado , Dictionary learning algorithms for sparse representation. Neural Comput. 15 , 349–396 (2003).12590811
2 I. Borg, P. Groenen, Modern Multidimensional Scaling: Theory and Applications. Springer Series in Statistics (2005).
3 K. Engan, S. O. Aase, J. H. Husoy, “Method of optimal directions for frame design” in 1999 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings. ICASSP99 (Cat. No.99CH36258) (1999), vol. 5, pp. 2443–2446.
4 L. Borcea, G. Papanicolaou, C. Tsogka, Interferometric array imaging in clutter. Inverse Probl. 21 , 1419–1460 (2005).
5 J. Aspnes , A theory of network localization, mobile computing. IEEE Trans. 5 , 1663–1678 (2007).
6 R. G. Stansfield, Statistical theory of D.F. fixing. J. Instit. Elect. Eng. 94 , 762–770 (1947).
7 A. Yeredor, E. Angel, Joint TDOA and FDOA estimation: A conditional bound and its use for optimally weighted localization. IEEE Trans. Signal Process. 59 , 1612–1623 (2011).
8 P. Wu , Time difference of arrival (TDoA) localization combining weighted least squares and firefly algorithm. Sensors 19 , 2554 (2019).31167498
9 Y. Shang, W. Ruml, Y. Zhang, M. P. J. Fromherz, “Localization from mere connectivity” in MobiHoc 2003: Proceedings of the 4th ACM International Symposium on Mobile Ad Hoc Networking & Computing (ACM, New York, NY, USA, 2003), pp. 201–212.
10 S. Oh, A. Montanari, A. Karbasi, “Sensor network localization from local connectivity: Performance analysis for the MDS-MAP algorithm” in 2010 IEEE Information Theory Workshop on Information Theory (ITW 2010, Cairo), (Cairo, Egypt, 2010), pp. 1–5.
11 J. Garnier, G. Papanicolaou, Passive Imaging with Ambient Noise (Cambridge University Press, United States, New York, NY, 2016).
12 A. Bakulin, R. Calvert, The virtual source method: Theory and case study. Geophysics 71 , SI139–SI150 (2006).
13 A. Agarwal, A. Anandkumar, P. Jain, P. Netrapalli, Learning sparsely used overcomplete dictionaries via alternating minimization. SIAM J. Optim. 26 , 2775–2799 (2016).
14 M. Fink , Time-reversed acoustics. Rep. Prog. Phys. 63 , 1933–1995 (2000).
15 L. Borcea, G. Papanicolaou, C. Tsogka, J. Berryman, Imaging and time reversal in random media. Inverse Probl. 18 , 1247–1279 (2002).
16 M. E. M. Aharon, M. Elad, A. Bruckstein, K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation. IEEE Trans. Sig. Process. 54 , 4311–4322 (2006).
17 A. Agarwal, A. Anandkumar, P. Netrapalli, A clustering approach to learning sparsely used overcomplete dictionaries. IEEE Trans. Inf. Theory 63 , 575–592 (2017).
18 A. Novikov, S. White, “Spectral subspace dictionary learning” in Proceedings of Machine Learning Research, 34th International Conference on Algorithmic Learning Theory, S. Agrawal, F. Orabona, Eds. (2023), pp. 1–36.
19 G. Davis, S. Mallat, M. Avellaneda, Adaptive Greedy approximations. J. Const. Approx. 13 , 57–98 (1997).
20 M. R. Osborne, B. Presnell, B. A. Turlach, On the lasso and its dual. J. Comput. Graph. Stat. 9 , 319–337 (2000).
21 I. Daubechies, M. Defrise, C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Comm. Pure Appl. Math. 57 , 1413–1457 (2004).
22 A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Img. Sci. 2 , 183–202 (2009).
23 M. Moscoso, A. Novikov, G. Papanicolaou, L. Ryzhik, A differential equations approach to l1-minimization with applications to array imaging. Inverse Probl. 28 , 105001 (2012).
24 L. Borcea, J. Garnier, G. Papanicolaou, C. Tsogka, Enhanced statistical stability in coherent interferometric imaging. Inverse Probl. 27 , 085004 (2011).
25 M. Moscoso, A. Novikov, G. Papanicolaou, C. Tsogka, Multifrequency interferometric imaging with intensity-only measurements. SIAM J. Img. Sci. 10 , 1005–1032 (2017).
26 S. M. Rytov, Y. A. Kravtsov, V. I. Tatarskii, “Principles of statistical radiophysics” in 4. Wave Propagation Through Random Media (Springer Verlag, Berlin, 1989).
27 M. Born, E. Wolf, Principles of Optics (Academic Press, 1970).
