Rapid evaluation of the spectral signal detection threshold and Stieltjes transform

TL;DR

This paper presents fast algorithms for evaluating spectral signal detection thresholds and the Stieltjes transform in large data matrices with heteroscedastic noise, aiding spectral denoising in statistical applications.

Contribution

The authors develop novel, provably correct algorithms for computing spectral thresholds and Stieltjes transforms in complex noise models, improving speed and accuracy over existing methods.

Findings

Algorithms are faster and more accurate than previous methods.

Validated performance through numerical experiments.

Applicable to models with separable variance profiles.

Abstract

Accurate detection of signal components is a frequently-encountered challenge in statistical applications with low signal-to-noise ratio. This problem is particularly challenging in settings with heteroscedastic noise. In certain signal-plus-noise models of data, such as the classical spiked covariance model and its variants, there are closed formulas for the spectral signal detection threshold (the largest sample eigenvalue attributable solely to noise) for isotropic noise in the limit of infinitely large data matrices. However, more general noise models currently lack provably fast and accurate methods for numerically evaluating the threshold. In this work, we introduce a rapid algorithm for evaluating the spectral signal detection threshold in the limit of infinitely large data matrices. We consider noise matrices with a separable variance profile (whose variance matrix is rank one), as these arise often in applications. The solution is based on nested applications of Newton's method. We also devise a new algorithm for evaluating the Stieltjes transform of the spectral distribution at real values exceeding the threshold. The Stieltjes transform on this domain is known to be a key quantity in parameter estimation for spectral denoising methods. The correctness of both algorithms is proven from a detailed analysis of the master equations characterizing the Stieltjes transform, and their performance is demonstrated in numerical experiments.

Figures

Tables

Table 1. Table 1 : Error after each iterate in evaluation of e λ ∗ superscript subscript 𝑒 𝜆 e_{\lambda}^{*} .

$m$	$T_{\lambda }\left(e_m\right)$	err($e_m,e_{\lambda }^{\ast }$)
1	2.08e-01	4.26e-02
2	4.34e-02	1.11e-02
3	2.55e-03	6.90e-04
4	9.66e-06	2.63e-06
5	1.40e-10	3.80e-11
6	-2.75e-16	-1.23e-16

Table 2. Table 2 : Error after each iterate in evaluation of t ​ ( λ ) 𝑡 𝜆 t(\lambda) .

$m$	$R_{\lambda }\left(e_m\right)$	err($e_m,t\left(\lambda \right)$)
1	-3.77e-01	3.71e-02
2	-6.89e-02	8.34e-03
3	-3.40e-03	4.34e-04
4	-9.24e-06	1.18e-06
5	-6.84e-11	8.73e-12
6	-5.55e-15	7.85e-16

Table 3. Table 3 : Error after each iterate in evaluation of λ ∗ superscript 𝜆 \lambda^{*} .

$m$	$Q\left({\lambda }_m\right)$	err(${\lambda }_m,{\lambda }^{\ast }$)
1	1.25e+00	-3.10e-01
2	3.10e-01	-1.03e-01
3	3.12e-02	-1.15e-02
4	3.92e-04	-1.46e-04
5	6.34e-08	-2.37e-08
6	2.00e-15	-2.11e-16

Table 4. Table 4 : Error after each iterate in evaluation of e ​ ( λ ) 𝑒 𝜆 e(\lambda) .

$m$	$F_{\lambda }(e_m$)	err($e_m,e\left(\lambda \right)$)
1	4.71e-01	-1.00e+00
2	8.86e-03	-1.93e-02
3	6.42e-06	-1.40e-05
4	3.43e-12	-7.50e-12
5	4.44e-16	-1.10e-15

Table 5. Table 5 : Timings in seconds for evaluating λ ∗ superscript 𝜆 \lambda^{*} and s ​ ( λ ) 𝑠 𝜆 s(\lambda) at 100 values of λ 𝜆 \lambda .

${\log }_2\left(n\right)$	Timing, ${\lambda }_{\ast }$	Timing, $s\left(\lambda \right)$
18	4.30e-01	1.93e-01
19	8.67e-01	4.27e-01
20	1.89e+00	1.66e+00
21	3.95e+00	4.20e+00
22	8.41e+00	1.02e+01
23	1.70e+01	2.09e+01
24	3.43e+01	4.19e+01

Table 6. Table 6 : Errors and bias in estimating λ ∗ superscript 𝜆 \lambda^{*} .

${\log }_2\left(k\right)$	Error	Bias
5	8.73e-02	7.72e-02
6	5.41e-02	4.75e-02
7	3.39e-02	2.93e-02
8	2.12e-02	1.82e-02
9	1.32e-02	1.13e-02
10	8.20e-03	6.97e-03
11	5.12e-03	4.32e-03

Table 7. Table 7 : Errors in estimating λ 𝜆 \lambda , c 𝑐 c , and c ¯ ¯ 𝑐 \underline{c} .

${\log }_2\left(k\right)$	sing. val.	left cos.	right cos.
5	7.68e-02	1.80e-02	2.26e-02
6	5.43e-02	1.27e-02	1.59e-02
7	3.85e-02	9.10e-03	1.13e-02
8	2.74e-02	6.38e-03	7.98e-03
9	1.92e-02	4.53e-03	5.61e-03
10	1.36e-02	3.20e-03	3.99e-03
11	9.58e-03	2.25e-03	2.81e-03

Table 8. Table 8 : Timings in seconds of Algorithm 3 and the solution via minimizing z ​ ( v ) 𝑧 𝑣 z(v) .

${\log }_2\left(p\right)$	Alg. 3	Min. $z\left(v\right)$
18	2.67e-01	2.38e-02
19	5.29e-01	4.63e-02
20	1.14e+00	1.30e-01
21	2.43e+00	3.20e-01
22	5.37e+00	6.94e-01
23	1.09e+01	1.38e+00
24	2.18e+01	2.73e+00

Equations

$${\lambda }_j=\{\begin{array}{ll}\left({\ell }_j+1\right)\left(1+\frac{\gamma }{{\ell }_j}\right),& \text{ if }{\ell }_j>\sqrt{\gamma }\\ (1+\sqrt{\gamma }{)}^2,& \text{ if }{\ell }_j\le \sqrt{\gamma }.\end{array}.$$

$$c_{jm}=\{\begin{array}{ll}\sqrt{\frac{{\ell }_j^2-\gamma }{{\ell }^2+\gamma {\ell }_j}}& \text{ if }j=m\text{ and }{\ell }_j>\sqrt{\gamma }\\ 0,& \text{ otherwise }\end{array}.$$

$$d{\mu }_k(t)=\frac{1}{k}\sum _{i=1}^k{\delta }_{{\lambda }_i}(t).$$

$$\gamma =\underset{k\to \infty }{\lim }\frac{k}{l(k)}$$

$${\lambda }^{\ast }=\mathrm{argmax}\{\lambda >0:\mu ([\lambda ,\infty ))>0\}.$$

$$s(\lambda )={\int }_{\mathbb{R}}\frac{1}{t-\lambda }d\mu (t),$$

$$s^{\prime }(\lambda )={\int }_{\mathbb{R}}\frac{1}{(t-\lambda {)}^2}d\mu (t).$$

$$\underline{s}(\lambda )={\int }_{\mathbb{R}}\frac{1}{t-\lambda }d\underline{\mu }(t),$$

$${\underline{s}}^{\prime }(\lambda )={\int }_{\mathbb{R}}\frac{1}{(t-\lambda {)}^2}d\underline{\mu }(t).$$

$$d\nu =\sum _{i=1}^p{\omega }_i{\delta }_{a_i}$$

$$d\underline{\nu }=\sum _{j=1}^n{\pi }_j{\delta }_{b_j}.$$

$$\underline{s}(\lambda )=\gamma s(\lambda )+\frac{\gamma -1}{\lambda }$$

$$\displaystyle s(\lambda)=\frac{1}{\gamma}\underline{s}(\lambda)+\bigg{(}\frac{1}{\gamma}-1\bigg{)}\frac{1}{\lambda}.$$

$${\underline{s}}^{\prime }(\lambda )=\gamma s^{\prime }(\lambda )+\frac{1-\gamma }{{\lambda }^2}$$

$$\displaystyle s^{\prime}(\lambda)=\frac{1}{\gamma}\underline{s}^{\prime}(\lambda)+\bigg{(}1-\frac{1}{\gamma}\bigg{)}\frac{1}{\lambda^{2}}.$$

$$d\underline{\mu }(t)=\gamma d\mu (t)+(1-\gamma ){\delta }_0(t).$$

$$s(\lambda )={\int }_{\mathbb{R}}\frac{1}{aG(e(\lambda ))-\lambda }d\nu (a),$$

$$G(e)={\int }_{\mathbb{R}}\frac{b}{1+\gamma be}d\underline{\nu }(b)$$

$$e(\lambda )={\int }_{\mathbb{R}}\frac{a}{aG(e(\lambda ))-\lambda }d\nu (a).$$

$$d\nu =\sum _{i=1}^p{\omega }_i{\delta }_{a_i},$$

$$d\underline{\nu }=\sum _{j=1}^n{\pi }_j{\delta }_{b_j},$$

$$s(\lambda )=\sum _{i=1}^p\frac{{\omega }_i}{a_{i}G(e(\lambda ))-\lambda }$$

$$e(\lambda )=\sum _{i=1}^p\frac{a_i{\omega }_i}{a_{i}G(e(\lambda ))-\lambda }$$

$$G(e)=\sum _{j=1}^n\frac{b_j{\pi }_j}{1+\gamma b_{j}e}.$$

$$\mathbf{Y}=\mathbf{X}+\mathbf{N},$$

$$D(\lambda )=\lambda s(\lambda )\underline{s}(\lambda ),$$

$${\theta }_m^2=\underset{k,l\to \infty }{\lim }\frac{1}{D({\lambda }_m)}.$$

$$\underset{k,l\to \infty }{\lim }|⟨{\mathbf{u}}_m,{\hat{\mathbf{u}}}_m⟩{|}^2=\frac{s({\lambda }_m)D({\lambda }_m)}{D^{\prime }({\lambda }_m)}$$

$$\underset{k,l\to \infty }{\lim }|⟨{\mathbf{v}}_m,{\hat{\mathbf{v}}}_m⟩{|}^2=\frac{\underline{s}({\lambda }_m)D({\lambda }_m)}{D^{\prime }({\lambda }_m)}$$

$$Y_j=X_j+b_j^{1/2}{\mathbf{A}}^{1/2}{G}_j,1\le j\le n,$$

Full text

In this work, we introduce a rapid algorithm for evaluating the spectral signal detection threshold in the limit of infinitely large data matrices. We consider noise matrices with a separable variance profile (whose variance matrix is rank one), as these arise often in applications. The solution is based on nested applications of Newton’s method. We also devise a new algorithm for evaluating the Stieltjes transform of the spectral distribution at real values exceeding the threshold. The Stieltjes transform on this domain is known to be a key quantity in parameter estimation for spectral denoising methods. The correctness of both algorithms is proven from a detailed analysis of the master equations characterizing the Stieltjes transform, and their performance is demonstrated in numerical experiments.

1 Introduction

Random matrix theory is an increasingly popular tool for deriving methods in statistical signal processing applications. Recent applications include signal denoising [41], [18], [29], [33], [31], covariance estimation [16], [31], filter correction in signal and image processing [4], [14], multiple-input multiple-output (MIMO) wireless communication [8], [5], and factor analysis [15], [13], to name just a few.

Part of the appeal of random matrix theory is that for many classes of random matrices, certain random quantities of statistical significance converge to well-defined limits as the matrix size grows infinitely large. By way of introduction, we illustrate this phenemonon in the spiked covariance model, introduced by Johnstone in [24]. Here, the user observes $l$ iid random vectors $Y_{1},\dots,Y_{l}$ in $\mathbb{R}^{k}$ of the form $Y_{j}=X_{j}+\varepsilon_{j}$, where $X_{j}$ (the signal component) is a random vector with covariance $\mathbf{\Sigma}$ of rank $r\ll\min(k,l)$ and eigenvalues $\ell_{1}>\dots>\ell_{r}>0$, and $\varepsilon_{j}$ (the noise component) is a zero mean Gaussian vector with covariance $\mathbf{I}_{k}$. Typically, the user wishes to estimate quantities related to the signal covariance $\mathbf{\Sigma}$, such as its rank, eigenvectors, and eigenvalues, from the signal-plus-noise observations $Y_{1},\dots,Y_{l}$.

It is well-known that when $k$ and $l$ are both large, there emerge simple relationships between the sample covariance $\widehat{\mathbf{\Sigma}}=\frac{1}{l}\sum_{j=1}^{l}(Y_{j}-\overline{Y})(Y_{j}-\overline{Y})^{T}$ (where $\overline{Y}=\frac{1}{l}\sum_{j=1}^{l}Y_{j}$ denotes the sample mean) and the signal covariance $\mathbf{\Sigma}$. For example, the largest $r$ eigenvalues of $\widehat{\mathbf{\Sigma}}$ converge almost surely to the following limits as $k,l\to\infty$ and $k/l\to\gamma>0$:

[TABLE]

Furthermore, the inner products $\langle\mathbf{u}_{j},\hat{\mathbf{u}}_{m}\rangle$ between the top $r$ eigenvectors $\mathbf{u}_{1},\dots,\mathbf{u}_{r}$ of $\mathbf{\Sigma}$ and the top $r$ eigenvectors $\hat{\mathbf{u}}_{1},\dots,\hat{\mathbf{u}}_{r}$ of $\widehat{\mathbf{\Sigma}}$ converge almost surely to the following limits:

[TABLE]

Proofs of these relationships may be found in the seminal work [38], while extensions to more general noise models can be found in [3]. These relationships have been used to devise estimators of $\mathbf{\Sigma}$, as in the work [16], as well as to predict the signal vectors $X_{1},\dots,X_{l}$ themselves, as in [41], [18], [29], [30].

Of particular significance for the present work is the value $(1+\sqrt{\gamma})^{2}$. This is the limiting value of the operator norm of the sample covariance $\frac{1}{l}\sum_{j=1}^{l}\varepsilon_{j}\varepsilon_{j}^{T}$ of the noise component alone as $k,l\to\infty$ and $k/l\to\gamma$; see [19]. More generally, for any sequence of random noise matrices whose operator norms have an almost sure asymptotic limit as the dimensions grow to infinity, we refer to that limiting value as the spectral signal detection threshold, or SSDT for short. As is apparent for white noise from the formula (1), observed eigenvalues exceeding the SSDT may be attributed to the signal component, and used to extract information about the signal. Eigenvalues below the SSDT are not as useful for recovering information on the signal (though see Remark 4 below). In typical applications of the spiked covariance model and its generalizations to signal denoising, only the eigenvalues of $\widehat{\mathbf{\Sigma}}$ exceeding the SSDT are used for denoising [18], [16], [30], [33], [14], [31], [4]. For this reason, evaluation of the SSDT is a critical task in statistical signal processing.

While the SSDT and related quantities, such as the relations between the eigenvalues of $\mathbf{\Sigma}$ and $\widehat{\mathbf{\Sigma}}$, have simple formulas in the case of white noise, more general noise models pose greater challenges. In this paper, we study a fairly broad family of noise matrices, namely those with a separable variance profile: every entry of the noise matrix has a potentially different variance, but the matrix of all the variances is rank one. This type of random matrix arises naturally in a number of applications, such as signal denoising with variable-strength heteroscedastic noise and the Kronecker model of multiple-input multiple-output transmission. We show that even when there are not closed formulas for the quantities of interest in this model, they may nevertheless be evaluated rapidly and to machine precision via provably correct algorithms.

1.1 This paper’s contributions

This paper introduces fast, scalable, and numerically stable algorithms for two related problems arising in statistical signal processing applications. The first is to evaluate the spectral signal detection threshold (SSDT) for noise matrices with a separable variance profile. As explained above, the SSDT is the asymptotic value of the operator norm of a random matrix representing the noise component in a signal-plus-noise observation model as the dimensions of the matrix grow to infinity.

The second problem we address is evaluating the Stieltjes transform (also known as the Cauchy-Stieltjes transform) of a certain probability measure associated with the noise matrix, known as the limiting spectral distribution (LSD). The Stieltjes transform is a function of central interest in random matrix theory and its applications; we recall its precise definition in Section 2.1. As we will review in Section 2.3, the Stieltjes transform, evaluated at real values exceeding the SSDT, may be used to estimate certain model parameters and perform signal denoising. Like the SSDT itself, in isotropic noise the Stieltjes transform has a closed form, while more general variance profiles pose greater challenges which we address in the current paper.

The algorithms we present scale linearly with the number of problem parameters, and are provably accurate to machine precision. The solution and analysis of each problem makes use of the master equations characterizing the Stieltjes transform, which will be reviewed in Section 2.2. Through a detailed analysis, presented in Section 3, we prove that the desired values may be computed by the use of Newton’s root-finding algorithm. For the SSDT, several intermediate quantities used in each iteration of Newton’s method are themselves computable by Newton’s method, as we will show. Consequently, all parameters we compute are either available analytically, or can be provably computed to full precision by a rapidly converging iterative scheme. For matrices of size $p$-by-$n$, the asymptotic costs of the algorithms scale like $O(p+n)$.

1.2 Relation to prior work

The signal-plus-noise matrix models of the kind we consider have been previously studied in the statistical literature, in works such as [24], [3], [33], [14], [31], [18], [38], [30], [16], [22], [20], [49], [21]. Estimating the number of signal terms in such principal components analysis and factor models is a well-studied problem in statistics and statistical signal processing applications [15], [13], [6], [23], [25], [37], [26]. Detection in the low SNR regime constitutes a distinct though conceptually relevant class of problems in signal processing [44], [46], [47], [48]. Noise matrices with separable variance profiles have also been studied in the context of the Kronecker product model for MIMO wireless communication [8], [5] and factor analysis with applications to economics [35]. We also note that such random matrices are a special case of “algebraic” random matrices, introduced in [40].

The solution to both of the problems we study rests on a known characterization of the Stieltjes transform of the LSD as the solution to a set of certain non-linear equations; these have appeared in [39] and [9]. We present a new, detailed analysis of these equations. As a consequence of our analysis, we show that evaluating the Stieltjes transform may be done by a straightforward application of Newton’s root-finding algorithm, whereas the SSDT is computable by several nested applications of Newton’s agorithm, appropriately initialized.

Previous works have considered the problem of evaluating the LSD, its Stieltjes transform, and the boundary of its support, for different classes of random noise matrices. The paper [12] proposes a scheme for evaluating the Stieltjes transform at complex values outside the support of the LSD; this method is based on a non-linear equation satisfied by the Stieltjes transform [42], [43], [32]. The method in [7] is devoted to extending the approach to mixture models. The papers [27], [28] are concerned with solving the inverse problem of recovering the population distribution from the observed spectrum; this problem is also taken up in [17]. The paper [15] contains a method for finding the boundary of a certain family of LSDs based on root-finding, although the use of Newton’s method is not employed or analyzed. To the best of our knowledge, the problem of evaluating the SSDT and Stieltjes transform for random matrices with a separable variance profile has not been addressed in prior work; and for special cases that have been studied, such as in [15], the prior work does not contain fast algorithms with the guarantees presented here.

1.3 Outline

The remainder of the paper is outlined as follows. In Section 2, we precisely state the problems we solve in this paper, and review the mathematical and numerical material that we will be using. In Section 3, we derive the core mathematical theory on which our algorithms rest, namely a detailed analysis of the master equations characterizing the Stieltjes transform of the LSD. In Section 4, we provide explicit descriptions of the numerical algorithms for finding the SSDT and evaluating the Stieltjes transform. In Section 5, we present the results of several illustrative numerical experiments demonstrating the performance of our algorithms.

2 Preliminaries

2.1 Setting and problem formulation

We suppose we have a $k$-by-$l$ random matrix of the form $\mathbf{N}=\mathbf{A}^{1/2}\mathbf{G}\mathbf{B}^{1/2}$, where $\mathbf{G}$ is a random matrix with iid entries of mean zero and variance $l^{-1}$, and $\mathbf{A}$ and $\mathbf{B}$ are positive-definite matrices of sizes $k$-by-$k$ and $l$-by-$l$, respectively. We define the empirical spectral distribution (ESD) $\mu_{k}$ to be the distribution of eigenvalues $\lambda_{1},\dots,\lambda_{k}$ of the matrix $\mathbf{N}\mathbf{N}^{T}$:

[TABLE]

As is typical in high-dimensional problems, we work in the asymptotic regime where $k$ and $l$ both grow to infinity. More precisely, we suppose that $l=l(k)$ grows with $k$, and that the limit

[TABLE]

is well-defined, positive, and finite. Furthermore, as $k$ and $l$ grow, we suppose that the eigenvalue spectra of $\mathbf{A}$ and $\mathbf{B}$ have well-defined asymptotic distributions $\nu$ and $\underline{\nu}$, respectively (in terms of weak convergence), with compact support. Under suitable conditions on $\mathbf{A}$ and $\mathbf{B}$, the sequence of measures $\mu_{k}$ will almost surely converge weakly to a compactly supported measure $\mu$, called the limiting spectral distribution (LSD). We also denote by $\underline{\mu}$ the limiting spectral distribution of the eigenvalues of $\mathbf{N}^{T}\mathbf{N}$.

We define the spectral signal detection threshold (SSDT) as follows:

[TABLE]

This is the right endpoint of the LSD $\mu$’s support, or equivalently the asymptotic operator norm of the matrix $\mathbf{N}\mathbf{N}^{T}$. As we will explain in Section 2.3, in a signal-plus-noise model $\mathbf{Y}=\mathbf{X}+\mathbf{N}$ where $\mathbf{X}$ is a low-rank signal matrix, eigenvalues exceeding $\lambda^{*}$ are attributable to the signal $\mathbf{X}$.

Next, we define the Stieltjes transform of $\mu$:

[TABLE]

which has derivative equal to

[TABLE]

Similarly, the Stieltjes transform of $\underline{\mu}$ is given by

[TABLE]

with derivative equal to

[TABLE]

We call $\underline{s}(\lambda)$ the associated Stieltjes transform of $\mu$. The functions $s$ and $\underline{s}$ are defined for all complex $\lambda$ outside the supports of $\mu$ and $\underline{\mu}$, respectively. However, as we will explain in Section 2.3, in the present work we are only interested in evaluating $s(\lambda)$, $\underline{s}(\lambda)$, and their derivatives for real $\lambda>\lambda^{*}$.

In this paper, we consider the setting where $\nu$ and $\underline{\nu}$ are discrete distributions of the form

[TABLE]

and

[TABLE]

Here, $a_{1},\dots,a_{p}$ and $b_{1},\dots,b_{n}$ are positive numbers, and the weights $\omega_{i}$ and $\pi_{j}$ satisfy $\sum_{i=1}^{p}\omega_{i}=\sum_{j=1}^{n}\pi_{j}=1$, $\omega_{i}>0$, and $\pi_{j}>0$. For example, if $k/2$ eigenvalues of $\mathbf{A}$ are 1, and the remaining $k/2$ are equal to 2, then $p=2$, $a_{1}=1$, $a_{2}=2$, and $\omega_{1}=\omega_{2}=1/2$.

With this background and notation, we can concisely state the two problems we solve in this paper:

Problem 1. Given $a_{1},\dots,a_{p}$ with corresponding weights $\omega_{1},\dots,\omega_{p}$; and $b_{1},\dots,b_{n}$ with corresponding weights $\pi_{1},\dots,\pi_{n}$; evaluate the SSDT $\lambda^{*}$.

Problem 2. Given $a_{1},\dots,a_{p}$ with corresponding weights $\omega_{1},\dots,\omega_{p}$; and $b_{1},\dots,b_{n}$ with corresponding weights $\pi_{1},\dots,\pi_{n}$; and a value $\lambda>\lambda^{*}$; evaluate $s(\lambda)$, $s^{\prime}(\lambda)$, $\underline{s}(\lambda)$, and $\underline{s}^{\prime}(\lambda)$.

Algorithms solving Problems 1 and 2 to machine precision, with asymptotic cost $O(p+n)$, are presented in Section 4.

Remark 1**.**

In the problem statements, the parameters $p$, $n$, $a_{1},\dots,a_{p}$, and $b_{1},\dots,b_{n}$ are user-supplied inputs. In many statistical applications, these parameters are assumed to be known, as in the works [15], [31], [12], [22], [20], [21].

Remark 2**.**

The assumption that $d\nu$ and $d\underline{\nu}$ are discrete measures is quite common in the literature on random matrix theory and its applications [15], [31], [12], [22], [20], [21]. It is likely that the methods in the present work can be generalized to more general measures, so long as certain linear functionals of the measures can be efficiently computed (e.g. by numerical quadrature). Pursuing this in detail lies outside the scope of the current work, however.

2.2 Properties of the Stieltjes transform

The Stieltjes and associated Stieltjes transforms of the LSD $\mu$ are defined by equations (6) and (8), respectively. We refer the reader to the standard references [1], [45] for a detailed treatment of Stieltjes transforms of probability measures, particularly random spectral measures.

Lemma 1**.**

The Stieltjes transforms $s(\lambda)$ and $\underline{s}(\lambda)$ satisfy the following relations:

[TABLE]

and

[TABLE]

Lemma 2**.**

The derivatives $s^{\prime}(\lambda)$ and $\underline{s}^{\prime}(\lambda)$ satisfy the following relations:

[TABLE]

and

[TABLE]

Lemma 2 follows from Lemma 1, which in turn follows immediately from the relation:

[TABLE]

Remark 3**.**

Lemmas 1 and 2 show that $\underline{s}(\lambda)$ and $\underline{s}^{\prime}(\lambda)$ may be easily evaluated once $s(\lambda)$ and $s^{\prime}(\lambda)$ are computed.

The next result is central to our subsequent analysis.

Theorem 1**.**

The Stieltjes transform $s(\lambda)$ for the LSD satisfies the following master equations:

[TABLE]

where

[TABLE]

and $e(\lambda)$ is a function that satisfies the equation

[TABLE]

For a proof of Theorem 1, see, for instance, the paper [39]. The paper [9] presents a slightly modified form of these equations, with a detailed analysis showing that $e(\lambda)$ is smooth for real $\lambda$ outside the support of $\mu$.

We assume that $\nu$ and $\underline{\nu}$ are discrete measures of the form

[TABLE]

and

[TABLE]

where $\sum_{i=1}^{p}\omega_{i}=\sum_{j=1}^{n}\pi_{j}=1$, $\omega_{i}>0$, and $\pi_{j}>0$. The master equations therefore become:

[TABLE]

where $e(\lambda)$ is a function satisfying

[TABLE]

and the function $G$ is defined by:

[TABLE]

We note that in the case where $\mathbf{B}=\mathbf{I}_{n}$ (the singly-weighted case), the master equations take on a simpler form; see, for instance, [43], [42].

2.3 Spiked models and the $D$-transform

In a spiked matrix model, we observe a signal-plus-noise matrix $\mathbf{Y}$ of the form

[TABLE]

where $\mathbf{X}=\sum_{m=1}^{r}\theta_{m}\mathbf{u}_{m}\mathbf{v}_{m}^{T}$ is a rank $r\ll\min\{k,l\}$ signal matrix, and $\mathbf{N}$ is a noise matrix. The $D$-transform is defined as follows:

[TABLE]

where $s(\lambda)$ and $\underline{s}(\lambda)$ are, respectively, the Stieltjes transform and associated Stieltjes transform of the LSD of $\mathbf{N}\mathbf{N}^{T}$. The $D$-transform was introduced in [3], though as a function of $\sqrt{\lambda}$ rather than $\lambda$.

Denoting by $\lambda_{1},\dots,\lambda_{r}$ the top $r$ eigenvalues of $\mathbf{Y}\mathbf{Y}^{T}$, and $\hat{\mathbf{u}}_{1},\dots,\hat{\mathbf{u}}_{r}$, $\hat{\mathbf{v}}_{1},\dots,\hat{\mathbf{v}}_{r}$ the corresponding left and right singular vectors of $\mathbf{Y}$, respectively, it is shown in [3] that the $D$-transform defines a mapping between $\lambda_{1},\dots,\lambda_{r}$ and the eigenvalues $\theta_{1}^{2},\dots,\theta_{r}^{2}$ of the signal matrix $\mathbf{X}\mathbf{X}^{T}$, which holds almost surely in the limit $k,l\to\infty$:

[TABLE]

This is satisfied for sufficiently large eigenvalues $\theta_{m}^{2}$ of $\mathbf{X}\mathbf{X}^{T}$, namely those for which $\theta_{m}^{2}>1/D(\lambda^{*})$. Phrased differently, an observed eigenvalue $\lambda_{m}$ of $\mathbf{Y}\mathbf{Y}^{T}$ satisfying $\lambda_{m}>\lambda^{*}$ may be attributed to the presence of signal.

Furthermore, the asymptotic cosines $\langle\mathbf{u}_{m},\hat{\mathbf{u}}_{m}\rangle$ and $\langle\mathbf{v}_{m},\hat{\mathbf{v}}_{m}\rangle$, $1\leq m\leq r$, can also be evaluated using the Stieltjes transform. It is shown that, almost surely,

[TABLE]

and

[TABLE]

Consequently, evaluating $s(\lambda_{m})$ and $s^{\prime}(\lambda_{m})$ provides a method for estimating the angles between the population singular vectors $\mathbf{u}_{m}$, $\mathbf{v}_{m}$ and the observed singular vectors $\hat{\mathbf{u}}_{m}$, $\hat{\mathbf{v}}_{m}$. These relationships have been employed to derive optimal methods of singular value shrinkage [33].

Remark 4**.**

While eigenvalues of $\mathbf{Y}\mathbf{Y}^{T}$ exceeding the SSDT $\lambda^{*}$ indicate the presence of signal, the converse statement – if all eigenvalues are less than $\lambda^{*}$, then there is no signal – is subtler. For example, the work [36] studies the detection of signals in white noise and shows that the joint distribution of all eigenvalues of $\mathbf{Y}\mathbf{Y}^{T}$ changes even when the signal eigenvalues do not exceed the SSDT. However, the resulting statistical test does not accurately distinguish between the presence or absence of signal with overwhelming probability. Regardless, in tasks such as covariance and signal estimation in the spiked model, only those eigenvalues exceeding the SSDT are typically used for covariance estimation and signal denoising [18], [16], [30], [33], [14], [31], [4].

Remark 5**.**

Separable variance profiles arise naturally when the columns of $\mathbf{Y}$ are independent random vectors in $\mathbb{R}^{p}$ of the form

[TABLE]

where the $X_{j}$ are signal vectors constrained to an $r$-dimensional subspace of $\mathbb{R}^{p}$; $\mathbf{A}$ is a positive definite matrix; and $b_{1},\dots,b_{n}$ are specified weights. In this model, each signal vector $X_{j}$ is observed in the presence of heteroscedastic noise (with covariance $\mathbf{A}$) of variable strength $b_{j}$. The noise matrix alone is distributed like $\mathbf{N}=\mathbf{A}^{1/2}\mathbf{G}\mathbf{B}^{1/2}$, where $\mathbf{B}=\text{diag}(b_{1},\dots,b_{n})$. Models like these are considered in, for example, [30], [11]. When $\mathbf{A}=\mathbf{I}_{p}$, this model is considered in [22], [20], [21]; when $\mathbf{B}=\mathbf{I}_{n}$, this model is considered in [49], [31].

Proposition 1**.**

For $\lambda>\lambda^{*}$, $D(\lambda)$ is decreasing and convex.

Proof.

Write $D(\lambda)=\varphi(\lambda)\underline{\varphi}(\lambda)$, where $\varphi(\lambda)=\sqrt{\lambda}s(\lambda)$ and $\underline{\varphi}(\lambda)=\sqrt{\lambda}\underline{s}(\lambda)$. Using $D^{\prime}(\lambda)=\varphi(\lambda)\underline{\varphi}^{\prime}(\lambda)+\varphi^{\prime}(\lambda)\underline{\varphi}(\lambda)$ and $D^{\prime\prime}(\lambda)=\varphi(\lambda)\underline{\varphi}^{\prime\prime}(\lambda)+\varphi^{\prime\prime}(\lambda)\underline{\varphi}(\lambda)+2\varphi^{\prime}(\lambda)\underline{\varphi}^{\prime}(\lambda)$, and $\varphi(\lambda)<0$ and $\underline{\varphi}(\lambda)<0$ for $\lambda>\lambda^{*}$, it is enough to show that $\varphi(\lambda)$ $\underline{\varphi}(\lambda)$ are increasing and concave. But this follows immediately from the identities

[TABLE]

and

[TABLE]

and the definitions (6) and (8) of $s(\lambda)$ and $\underline{s}(\lambda)$. ∎

2.4 Newton’s root-finding algorithm

Newton’s method is a classical technique for finding the roots of a function of one real variable. We briefly review the method here; the reader may consult any standard reference on optimization or numerical analysis, such as [34], [10], for additional details. We are given a smooth function $f(x)$, where $x\in[a,b]$, and we suppose $f(a)<0$ and $f(b)>0$. We also suppose that $f^{\prime}(x)>0$ and $f^{\prime\prime}(x)<0$ for all $x\in(a,b)$; that is, $f$ is a strictly increasing and concave. Our goal is to compute $x^{*}$, the unique root of $f$ in $(a,b)$.

To find $x^{*}$, Newton’s method initializes $a<x_{0}<x^{*}$, and defines a sequence of updates recursively as follows: given an estimate $x_{k}$, the next value $x_{k+1}$ is defined by

[TABLE]

Geometrically, $x_{k+1}$ is the root of the line tangent to the graph of $f$ at the point $(x_{k},f(x_{k}))$. Because $f$ is concave, it is easy to see that $x_{k}<x_{k+1}\leq x^{*}$.

Because each $x_{k}$ is obtained from a linear approximation to $f$, the errors decay quadratically in the vicinity of $x^{*}$; that is

[TABLE]

when $|x_{k}-x^{*}|$ is sufficiently small.

Remark 6**.**

Quadratic convergence is what makes Newton’s method an especially attractive algorithm when it is applicable – that is, when both $f(x)$ and $f^{\prime}(x)$ may be accurately evaluated, and $f$ exhibits the conditions specified above. In practical terms, quadratic convergence means that the number of accurately computed digits of $x^{*}$ approximately doubles after each iteration of the algorithm, until machine precision is reached.

3 Mathematical apparatus

In this section, we analyze the master equations (22) – (24). Our results will provide the necessary tools to devise algorithms for computing the SSDT $\lambda^{*}$ and evaluating $s(\lambda)$ and $s^{\prime}(\lambda)$ for $\lambda>\lambda^{*}$. We define

[TABLE]

and

[TABLE]

We define the function $F(\lambda,e)$ by:

[TABLE]

Then for each $\lambda>\lambda^{*}$, $e(\lambda)$ satisfies $F(\lambda,e(\lambda))=0$. When we treat $\lambda$ as a fixed parameter and $e$ as a variable, we will use the notation $F_{\lambda}(e)=F(\lambda,e)$.

3.1 Range and monotonicity of $e(\lambda)$ when $\lambda>\lambda^{*}$

We first state a result on the range of $e(\lambda)$ for $\lambda>\lambda^{*}$. The function $G(e)$ approaches [math] as $e\to\infty$, and grows to $+\infty$ as $e\to(-1/\gamma b^{*})^{+}$, and is strictly decreasing on the interval

[TABLE]

For any $\lambda>0$, we define the interval $I_{\lambda}$ by

[TABLE]

Proposition 2**.**

When $\lambda>\lambda^{*}$, $e(\lambda)$ is contained in the interval $I_{\lambda}\cap(-\infty,0)$.

We will develop the proof in several steps.

Lemma 3**.**

Let $\epsilon>0$. Then the function $G(e(\lambda))$ is bounded for all $\lambda>\lambda^{*}+\epsilon$.

Proof.

We show that the range of $e(\lambda)$ cannot approach any of the singularities of $G$, which lie at the values $-1/(\gamma b_{j})$. Define

[TABLE]

Then from (22)

[TABLE]

and consequently

[TABLE]

Since $\lambda s(\lambda)$ is bounded for $\lambda>\lambda^{*}+\epsilon$, this tells us that $G(e(\lambda))$ must stay bounded so long as $e(\lambda)$ is bounded away from [math]; in particular, $e(\lambda)$ cannot be made arbitrarily close to any of the singularities $-1/(\gamma b_{j})$. ∎

Corollary 1**.**

$e(\lambda)\to 0^{-}$* as $\lambda\to\infty$.*

Proof.

Because $G(e(\lambda))$ is bounded for large $\lambda$, the result follows from (23). ∎

Corollary 2**.**

For any $\lambda>\lambda^{*}$, $e(\lambda)\in J$; that is,

[TABLE]

Proof.

This follows from the continuity of $e(\lambda)$ for $\lambda>\lambda^{*}$, and the facts that it never passes through $-1/(\gamma b_{j})$ and approaches 0 at large $\lambda$. ∎

Corollary 3**.**

For all $\lambda>\lambda^{*}$, $G(e(\lambda))>0$ and $e(\lambda)<0$.

Proof.

The positivity of $G(e(\lambda))$ follows immediately from Corollary 2. The negativity of $e(\lambda)$ then follows from (42), and the fact that $1+\lambda s(\lambda)<0$. ∎

Lemma 4**.**

For all $\lambda>\lambda^{*}$,

[TABLE]

Proof.

We have $G(e(\lambda))\neq\lambda/a_{i}$, from (22). Suppose for some $\lambda>\lambda^{*}$, we had

[TABLE]

The inequality must then remain true for all sufficiently large $\lambda$, since $G(e(\lambda))/\lambda$ is continuous and does not pass through $1/a_{i}$. However, the left side converges to [math] as $\lambda\to\infty$, since $G(e(\lambda))$ is bounded for large $\lambda$; a contradiction. This completes the proof. ∎

We have shown that for all $\lambda>\lambda^{*}$, $e(\lambda)$ lies in the interval defined by the inequalities

[TABLE]

This completes the proof of Proposition 2.

Next we prove that $e(\lambda)$ is increasing:

Proposition 3**.**

The function $e(\lambda)$ is increasing for $\lambda>\lambda^{*}$.

Proof.

We have:

[TABLE]

Since $F(\lambda,e(\lambda))=0$, we have:

[TABLE]

and so

[TABLE]

Consequently, $e^{\prime}(\lambda)$ can never be [math]. Furthermore,

[TABLE]

which converges to 1 as $\lambda\to\infty$ (note that $e(\lambda)$ stays bounded away from singularities of $G(e)$ and also $G^{\prime}(e)$, which have the same singularities). So $e^{\prime}(\lambda)>0$ for sufficiently large $\lambda$, and hence, since it is continuous and cannot pass through 0, $e^{\prime}(\lambda)>0$ for all $\lambda>\lambda^{*}$. ∎

3.2 Behavior of $F_{\lambda}(e)$

In this section we characterize the behavior of $F_{\lambda}(e)=F(\lambda,e)$ (viewed as a function of $e$) on the interval $I_{\lambda}$. Specifically, we show the following. For any $\lambda>0$, $F_{\lambda}(e)$ is a strictly convex function that approaches $+\infty$ as $e$ approaches either end of $I_{\lambda}$. Furthermore, when $\lambda>\lambda^{*}$, the minimum value of $F_{\lambda}(e)$ is less than zero, and $F_{\lambda}(0)>0$; consequently, there are exactly two roots of $F_{\lambda}(e)$, both contained in the interval $I_{\lambda}\cap(-\infty,0)$. We show that $e(\lambda)$ is always equal to the largest root, namely the one at which $F_{\lambda}^{\prime}(e)>0$.

Lemma 5**.**

For $\lambda>0$, the function $F_{\lambda}(e)$ is strictly convex on $I_{\lambda}$; that is, $(\partial_{ee}^{2}F)(\lambda,e)>0$.

Proof.

We have:

[TABLE]

and

[TABLE]

Now whenever $e\in I_{\lambda}$, we have $e>-1/\gamma b^{*}\geq-1/\gamma b_{j}$ for all $1\leq j\leq n$, and $G(e)<\lambda/a^{*}\leq\lambda/a_{i}$ for all $1\leq i\leq p$; consequently,

[TABLE]

and

[TABLE]

Consequently, $(\partial_{ee}^{2}F)(\lambda,e)>0$ for all $e\in I_{\lambda}$, i.e. the function $F_{\lambda}(e)$ is convex. ∎

Lemma 6**.**

The function $F_{\lambda}(e)$ diverges to $+\infty$ as $e\to+\infty$ and as $e$ approaches the left endpoint of $I_{\lambda}$ from the right.

Proof.

This is immediate from the definition of $F(\lambda,e)$. ∎

Lemma 7**.**

For each $\lambda>\lambda^{*}$, $F(\lambda,0)>0$.

Proof.

We have:

[TABLE]

Since $G(e(\lambda))<\lambda/a_{i}$ and $G(e)$ is decreasing on $I_{\lambda}$, and $e(\lambda)<0$, we also have $G(0)<G(e(\lambda))<\lambda/a_{i}$; consequently, $F(\lambda,0)>0$. ∎

Proposition 4**.**

For $\lambda>\lambda^{*}$, $(\partial_{e}F)(\lambda,e(\lambda))>0$.

Proof.

This follows from (49) and $e^{\prime}(\lambda)>0$. ∎

We have shown that $F_{\lambda}(e)$ has two roots in the interval $I_{\lambda}\cap(-\infty,0)$ whenever $\lambda>\lambda^{*}$. Proposition 4 identifies $e(\lambda)$ as the root that is closest to zero, or equivalently, the root where the derivative of $F_{\lambda}$ is positive. This characterization will be used in Section 4.2 to devise the algorithm for computing $e(\lambda)$, and consequently $s(\lambda)$.

3.3 The minimum of $F_{\lambda}(e)$

We will let $t(\lambda)$ denote the minimum of $F_{\lambda}(e)$ on $J_{\lambda}$; that is, $t(\lambda)$ is the unique value on $J_{\lambda}$ that satisfies

[TABLE]

We define the function $Q(\lambda)$ for $\lambda>0$ by:

[TABLE]

For any $\lambda>0$, we define the function $R_{\lambda}(e)$ for $e\in I_{\lambda}$ by:

[TABLE]

Note that by definition, $R_{\lambda}(t(\lambda))=0$ for all $\lambda=0$.

We will show that $Q$ is decreasing and convex and $R_{\lambda}$ is increasing and concave.

Lemma 8**.**

$Q(\lambda)$* is a decreasing function of $\lambda>0$.*

Proof.

The derivative of $Q$ may be computed as follows:

[TABLE]

which is negative. ∎

Proposition 5**.**

$Q(\lambda)$* is a convex function of $\lambda>0$.*

Proof.

We first compute the derivative of $t(\lambda)$. We have

[TABLE]

Differentiating with respect to $\lambda$, we find

[TABLE]

and so

[TABLE]

Now the second derivative of $Q$ is given by

[TABLE]

We will show that $Q^{\prime\prime}(\lambda)>0$ for all $\lambda$. In fact, we will show the stronger result that each summand is positive, or equivalently, recalling that $a^{*}=\max_{1\leq i\leq p}a_{i}$,

[TABLE]

To prove this, we observe that

[TABLE]

To show that this is less than 1, it is enough to show that

[TABLE]

or equivalently

[TABLE]

We will show that this inequality holds for any value of $a_{i}$ between [math] and $a^{*}$. If $a_{i}=a^{*}$, then the right side becomes

[TABLE]

and since the term

[TABLE]

is positive (because $G$ is convex, $t(\lambda)\in J_{\lambda}$, and $a^{*}G(e)<\lambda$ on $J_{\lambda}$), the desired inequality is satisfied.

On the other hand, if $a_{i}=0$, the right hand side becomes

[TABLE]

where the inequality holds since $t(\lambda)\in J_{\lambda}$, and $a^{*}G(e)<\lambda$ for all $e\in J_{\lambda}$. By the Cauchy-Schwarz inequality,

[TABLE]

or in other words,

[TABLE]

This gives the desired result. ∎

Proposition 6**.**

For all $\lambda>0$, $R_{\lambda}(e)$ is an increasing and concave function of $e\in I_{\lambda}$.

Proof.

The first derivative of $R_{\lambda}$ is:

[TABLE]

which as we’ve seen is positive. The second derivative of $R_{\lambda}$ is

[TABLE]

which is always negative, since $G^{(3)}(e)<0$, $G^{\prime\prime}(e)>0$, $G^{\prime}(e)<0$, and $a_{i}G(e)<\lambda$ for all $e\in I_{\lambda}$ and $1\leq i\leq p$. ∎

4 Algorithms

In this section, we describe the algorithms for computing the SSDT $\lambda^{*}$ and for evaluating the Stieltjes transform $s(\lambda)$ and its derivative $s^{\prime}(\lambda)$ at values $\lambda>\lambda^{*}$. By Lemmas 1 and 2, $\underline{s}(\lambda)$ and $\underline{s}^{\prime}(\lambda)$ may be easily found from $s(\lambda)$ and $s^{\prime}(\lambda)$.

4.1 Computation of the boundary $\lambda^{*}$

In this section, we derive an algorithm for the computation of $\lambda^{*}$. First, we observe that when $\lambda>\lambda^{*}$, then as we have shown there are real roots of $F_{\lambda}(e)$ to the left and to the right of $t(\lambda)$; in particular, $F(\lambda,t(\lambda))<0$. On the other hand, if $\lambda<\lambda^{*}$, the function $F_{\lambda}(e)$ cannot have a real root, since that would imply that the Stieltjes transform is real inside the support of $\mu$. Consequently, $F(\lambda,t(\lambda))>0$. It follows that the SSDT $\lambda^{*}$ is the value at which $F(\lambda,t(\lambda))=0$; that is, $\lambda^{*}$, is the unique root of $Q$ on $(0,\infty)$.

From Lemma 8 and Proposition 5, $Q(\lambda)$ is decreasing and convex. With an efficient procedure for evaluating $Q(\lambda)$ and $Q^{\prime}(\lambda)$, we can therefore use Newton’s algorithm to find its root if we initialize the algorithm to the left of the root. In Section 4.1.2, we detail how to evaluate $Q(\lambda)$ and $Q^{\prime}(\lambda)$. As a preliminary step, we will need to compute the left endpoint of $I_{\lambda}$; we do this in 4.1.1. The resulting algorithm for evaluating $\lambda^{*}$ is summarized in Algorithm 3.

4.1.1 Computation of the left endpoint of $I_{\lambda}$, $\lambda>0$

We recall the definition of the interval $I_{\lambda}$:

[TABLE]

where $J$ is the interval $J=\left\{e:e>-1/(\gamma b^{*})\right\}$. Let us denote by $e_{\lambda}^{*}$ the left endpoint of $I_{\lambda}$. Since $G(e)$ is a decreasing function of $e$ on $J$, $e_{\lambda}^{*}$ is the unique root of

[TABLE]

on $J$. Since $T_{\lambda}$ is a decreasing, convex function of $e$ on $J$, we may find its root using Newton’s algorithm, initialized to the left of the root. Such an initial value $e_{0}$ can be found by starting with any value $e$ in $J$, and performing bisection with $-1/(\gamma b^{*})$, the left endpoint of $J$, until we arrive at a value where $T_{\lambda}(e_{0})>0$. Newton’s algorithm is then performed with initial value $e_{0}$. We summarize the procedure in Algorithm 1.

4.1.2 Evaluation of $t(\lambda)$, $Q(\lambda)$, and $Q^{\prime}(\lambda)$

Next, we show how to evaluate the functions $t(\lambda)$, $Q(\lambda)$, and $Q^{\prime}(\lambda)$. $t(\lambda)$ is defined as the root of $R_{\lambda}(e)$ on $I_{\lambda}$. Since $R_{\lambda}(e)$ is an increasing, concave function, we may find its root using the Newton algorithm initialized to the left of the root, i.e. the region where $R_{\lambda}(e)<0$. Such an initial value may be found by taking a starting point $e$ to the right of $e_{\lambda}^{*}$, and performing bisection on $e$ and $e_{\lambda}^{*}$ until we arrive at a value $e_{0}$ with $R_{\lambda}(e_{0})<0$. We can then perform Newton’s algorithm on $R_{\lambda}$, initialized at $e_{0}$. The procedure is summarized in Algorithm 2.

4.2 Evaluation of $s(\lambda)$ and $s^{\prime}(\lambda)$, $\lambda>\lambda^{*}$

In this section, we present an algorithm for evaluating the Stieltjes transform $s(\lambda)$ of $\mu$, and its derivative $s^{\prime}(\lambda)$, when $\lambda>\lambda^{*}$. This immediately provides a method for evaluating $\underline{s}(\lambda)$ and $\underline{s}^{\prime}(\lambda)$, and the $D$-transform $D(\lambda)$ defined in Section 2.3.

As we showed in Section 3.2, the function $F_{\lambda}(e)$ is convex on $I_{\lambda}$ and has two roots, both of which are negative. The root closest to [math], i.e. the rightmost root, is $e(\lambda)$. Since $F_{\lambda}(0)>0$ and $F_{\lambda}(e)$ is convex, this tells us that Newton’s method, initialized at $e_{0}=0$, will converge to $e(\lambda)$. For brevity, we introduce the function $W(\lambda,e)$ defined by:

[TABLE]

With this notation, we have:

[TABLE]

and

[TABLE]

Note that from (48), we can evaluate $e^{\prime}(\lambda)$:

[TABLE]

The method for evaluating $s(\lambda)$ and $s^{\prime}(\lambda)$ is summarized in Algorithm 4.

5 Numerical results

In this section, we report on several experiments illustrating the behavior of the algorithms from Section 4. For the timings reported in Sections 5.2 and 5.5, we used an implementation written in MATLAB 2019b and run on a Dell Precision 5540 with 62.5 GB of RAM and an Intel Core i9 CPU. The MATLAB code is available at the following URL: https://github.com/wleeb/MPBoundary.

5.1 Convergence

Algorithms 1 – 4 are all versions of the Newton root-finding method. In this section, we illustrate the quadratic convergence of these methods, as predicted from the theory reviewed in Section 2.4. In each experiment, we used parameters $p=512$, $n=1024$, and $\gamma=1/2$. We generated the values $a_{1},\dots,a_{p}$ and $b_{1},\dots,b_{n}$ randomly from a $\text{Unif}(0,1)$ distribution, and assigned random probabilities $\omega_{i}$ and $\pi_{j}$ by drawing values from $\text{Unif}(0,1)$ and normalizing to sum to 1. For experiments in which a value $\lambda$ is specified, we also choose it at random.

In Tables 1 – 4, the first column displays the iteration number $m$, starting from the initial value and going until the root has been reached. The second column displays the function value at the $m^{th}$ iterate; the algorithm terminates when the function reaches machine precision $\epsilon$. We work in double precision, so $\epsilon\approx 10^{-16}$. The third column displays the relative error in the root itself, defined by:

[TABLE]

Table 1 shows the results for Algorithm 1, which computes the left endpoint $e_{\lambda}^{*}$ of $I_{\lambda}$. Table 2 shows the results for Algorithm 2, which evaluates $t(\lambda)$. Table 3 shows the results for Algorithm 3, which computes the boundary $\lambda^{*}$. Table 4 shows the results for Algorithm 4, which evaluates $e(\lambda)$.

Remark 7**.**

For each algorithm, we observe quadratic convergence close to the root, as expected. That is, on each iteration close to termination the number of correct digits roughly doubles, and the size of the objective roughly squares, until machine precision is reached.

5.2 Scalability

In the next experiment, we compute timings for the computation of $\lambda^{*}$ and the evaluation of $s(\lambda)$. For increasing values of $n$, we set $p=n/2$ ($\gamma=1/2$). We generated the values $a_{1},\dots,a_{p}$ and $b_{1},\dots,b_{n}$ randomly from a $\text{Unif}(1,2)$ distribution, and assigned them random probabilities $\omega_{i}$ and $\pi_{j}$ by drawing values from $\text{Unif}(0,1)$ and normalizing to sum to 1.

For each $n$, we then record the time in seconds required to compute $\lambda^{*}$, and the time in seconds required to compute $s(\lambda)$ on a grid of 100 equispaced values of $\lambda$ between $\lambda^{*}+1$ and $\lambda^{*}+10$. The reported timings are averaged over five runs of the experiment, and are displayed in Table 5. It is apparent that the running times scale approximately linearly with $n$, as we would expect. In addition to linearly scaling with $n$, the magnitudes of the timings are quite encouraging; for example, it takes only about 4 seconds to compute $\lambda^{*}$ when $n$ is over two million and $p$ is over one million.

5.3 Finite sample accuracy for $\lambda^{*}$

The master equations (22) – (24) from which we compute $\lambda^{*}$ are asymptotic and deterministic, holding almost surely in the limit as $k,l\to\infty$. In this experiment, we assess the finite sample accuracy of estimating the SSDT $\lambda^{*}$ from the operator norm of the random matrix $\mathbf{N}\mathbf{N}^{T}$. We consider how the estimate improves as $k$ and $l$ grow. We use a model where $p=n=2$, and both $\nu$ and $\underline{\nu}$ have point masses at $2$ and $3$, each with weight $1/2$; and set $\gamma=1/2$. We generate matrices of size $k$-by-$l$, where $k$ grows and $l=k/\gamma=2k$.

For each value of $k$, we draw such a $k$-by-$l$ random matrix $\mathbf{N}$ and compute the operator norm of $\mathbf{N}\mathbf{N}^{T}$. We compare this to the value of $\lambda^{*}$ computed using Algorithm 3. In Table 6, we show the mean absolute error, averaged over $M=40000$ runs for each value of $k$. More precisely, if $\hat{\lambda}_{j}^{*}$ is the operator norm of $\mathbf{N}\mathbf{N}^{T}$ from trial $j=1,\dots,M$, we record

[TABLE]

We also record the average bias, defined as

[TABLE]

In Figure 1, we plot the log error against $\log_{2}(k)$. The plot demonstrates a linear dependence. The slope of the line is estimated to be approximately $-0.68$.

Remark 8**.**

In [24], it is shown that for white noise the expected fluctations of the top eigenvalue of $\mathbf{N}\mathbf{N}^{T}$ are of the order $k^{-2/3}$, which implies the log-log plot would have a slope of $-2/3$. The observed slope of approximately $-0.68$ is a close match to this value.

Remark 9**.**

For all values of $k$, the bias is positive. In other words, the estimated values $\hat{\lambda}_{k}^{*}$ tend to underestimate $\lambda^{*}$.

5.4 Finite sample accuracy of spiked model parameters

In this experiment, we test the finite sample accuracy of parameter estimation in the spiked random matrix model. We consider $k$-by-$l$ random matrices of the form $\mathbf{Y}=\mathbf{X}+\mathbf{N}$, where $\mathbf{X}=\theta\mathbf{u}\mathbf{v}^{T}$ is a rank $1$ “signal” matrix with uniformly random singular vectors $\mathbf{u}$ and $\mathbf{v}$, and $\mathbf{N}=\mathbf{A}^{1/2}\mathbf{G}\mathbf{B}^{1/2}$ is a random Gaussian noise matrix with separable variance profile. We will denote by $\hat{\lambda}$

As we reviewed in Section 2.3, the top eigenvalue of $\mathbf{Y}\mathbf{Y}^{T}$ converges almost surely to a value $\lambda$ satisfying $\theta^{2}=1/D(\lambda)$, in the limit $k/l\to\gamma$. Furthermore, if $\hat{\mathbf{u}}$ and $\hat{\mathbf{v}}$ are the top singular vectors of $\mathbf{Y}$, then the absolute inner products $|\langle\mathbf{u},\hat{\mathbf{u}}\rangle|$ and $|\langle\mathbf{v},\hat{\mathbf{v}}\rangle|$ converge almost surely to $c\equiv s(\lambda_{m})D(\lambda_{m})/D^{\prime}(\lambda_{m})$ and $\underline{c}\equiv\underline{s}(\lambda_{m})D(\lambda_{m})/D^{\prime}(\lambda_{m})$, respectively. We remind the reader that we define the $D$-transform by $D(\lambda)=\lambda s(\lambda)\underline{s}(\lambda)$.

We test the accuracy of these formulas for finite and increasing values of $k$ and $l$. Using Algorithm 4 for evaluating $s(\lambda)$ and $s^{\prime}(\lambda)$, we can easily evaluate $D(\lambda)$. Proposition 1 shows that for a specified parameter $\theta$, Newton’s root-finding algorithm may be used to solve for the asymptotic $\lambda=1/D^{-1}(\theta^{2})$. The asymptotic values $c$ and $\underline{c}$ are then evaluated from their respective formulas.

We compare these asymptotic values to the top eigenvalue of $\mathbf{N}\mathbf{N}^{T}$ and the cosines between the singular vectors of $\mathbf{X}$ and $\mathbf{Y}$, for randomly generated data. We again use a model where $p=n=2$, and both $\nu$ and $\underline{\nu}$ have point masses at $2$ and $3$, each with weight $1/2$; and set $\gamma=1/2$. We generate matrices of size $k$-by-$l$, where $k$ grows and $l=k/\gamma=2k$. We generate a signal matrix $\mathbf{X}=\theta\mathbf{u}\mathbf{v}^{T}$ where $\mathbf{u}$ and $\mathbf{v}$ are uniformly random unit vectors in $\mathbb{R}^{k}$ and $\mathbb{R}^{l}$, respectively; and $\theta=\sqrt{1/D(\lambda^{*})+20}$, ensuring a detectable signal.

For each value of $k$, we draw such a $k$-by-$l$ random matrix $\mathbf{Y}=\mathbf{X}+\mathbf{N}$ and compute the operator norm of $\mathbf{Y}\mathbf{Y}^{T}$. We compare this to the asymptotic value of $\lambda$. We also compare the values $|\langle\mathbf{u},\hat{\mathbf{u}}\rangle|$ and $|\langle\mathbf{v},\hat{\mathbf{v}}\rangle|$ to $c$ and $\underline{c}$, respectively. In Table 7, we show the mean absolute errors of these estimates, averaged over $M=40000$ runs for each value of $k$. More precisely, if $\hat{\lambda}_{j}$ is the operator norm of $\mathbf{Y}\mathbf{Y}^{T}$ from trial $j=1,\dots,M$, and $\hat{c}_{j}$ and $\hat{\underline{c}}_{j}$ are the left and right cosines, respectively, we record the mean absolute errors:

[TABLE]

In Figure 2, we plot the log errors against $\log_{2}(k)$. The plots all demonstrate a linear dependence. The slope of each line is estimated to be approximately $-0.50$.

Remark 10**.**

In [3], it is shown that the expected fluctations of the top eigenvalue of $\mathbf{Y}\mathbf{Y}^{T}$ are of the order $k^{-1/2}$, which implies the log-log plot would have a slope of $-1/2$. In [2], it is shown that in the case of white noise the expected fluctations of the cosines are also of order $k^{-1/2}$, also leading to a slope of $-1/2$. Our observed slopes closely match these values.

5.5 One-sided weights

While Algorithm 3 computes the SSDT for an arbitrary separable variance profile, a special case of this problem that arises in certain applications is when $\mathbf{B}=\mathbf{I}_{n}$; that is, the variance profile has one-sided weights. In [15], it is shown that the SSDT may be evaluated by finding the minimizer $v^{*}$ of the function

[TABLE]

on the interval $(-1/a^{*},0)$, and then setting $\lambda_{*}=z(v^{*})$. The minimizer $v^{*}$ is the unique root of the monotonic function

[TABLE]

on the interval $(-1/a^{*},0)$.

When viable, this approach has obvious advantages over Algorithm 3, namely that it finds the root of a function which can be evaluated in closed form rather than by nested applications of Newton’s method. Properly applied, it should be substantially faster than Algorithm 3. However, the paper [15] does not analyze the behavior of the function (86), beyond observing that it is monotonic.

We compare the method from [15] for evaluating $\lambda_{*}$ based on the minimizer of $z(v)$ to Algorithm 3. To find the root of $z^{\prime}(v)$, we employ bisection until the error is approximately square root of machine epsilon, and then use Newton’s method to achieve full accuracy. In our experiment, we set $\gamma=1/2$, and for each value of $p$ we generate the $a_{i}$’s uniformly randomly on $[0,1]$ and take uniform weights $\omega_{i}=1/p$. For each value of $p$, we solve the problem using each method for $10$ runs, and average the timings of all the runs; the timings are presented in Table 8.

The results of the experiment are extremely encouraging. They suggest that for one-sided weights, working with the function (86) yields a substantially faster computation of the SSDT than the general master equations (22) – (24). However, a more detailed analysis of the function (86) must be carried out to justify such an approach and show that a fast algorithm is viable for general $a_{i}$ and $\omega_{i}$; this is beyond the scope of the current work.

6 Conclusion

We have introduced an algorithm for rapidly computing the spectral signal detection threshold $\lambda^{*}$ in signal-plus-noise random matrix models. We have considered a class of random noise matrices with separable variance structure, which arise often in applications. Our algorithm is based on an implicit characterization of $\lambda^{*}$ derived from the master equations for the Stieltjes transform $s(\lambda)$. Several nested applications of Newton’s method are applied to evaluate $\lambda^{*}$ and auxiliary parameters. Fast algorithms are also introduced to evaluate $s(\lambda)$ and $s^{\prime}(\lambda)$, the Stieltjes transform and its derivative, at real values $\lambda>\lambda^{*}$. We have demonstrated the rapid convergence of these methods and their linear scaling in numerical tests. We have also shown the effects of random fluctations for increasing values of $p$ and $n$.

Acknowledgements

I acknowledge support from NSF BIGDATA award IIS 1837992 and BSF award 2018230. I thank Edgar Dobriban for helpful discussions and for pointing out the method from [15]. I also thank the reviewers for their helpful comments on the manuscript.