# Robust estimation of the variance with STARRS

#### Daphné Giorgi, Antoine Godichon-Baggioni, Stéphane Robin, Laure Sansonnet

# 1 Framework

In all the sequel, we suppose that \(X\) is a random vector lying in \(\mathbb{R}^{p}\) with mean \(\mu\) and positive variance \(\Sigma\). Let us recall that the geometric median \(m\) of \(X\) is defined as the minimizer of the convex functional \(G : \mathbb{R}^p \to \mathbb{R}\), given for all \(h \in \mathbb{R}^p\) by Haldane (1948):

\[ G(h) := \mathbb{E} \left[ \| X - h \| - \| X \| \right], \]

where \(\| \cdot \|\) denotes the Euclidean norm.The median can be computed using one of two available algorithms, implemented respectively in the functions `WeiszfeldMedian` and `ASGMedian`. If the distribution of \(X\) is symmetric, then the mean \(\mu\) of \(X\) coincides with its median \(m\). Unfortunately, there is no direct robust alternative of the variance. Nevertheless one can focus on the Median Covariance Matrix, defined as the minimizer of the functional \(G_m : \mathcal{M}_p(\mathbb{R}) \to \mathbb{R}\), given for any square matrix \(M \in \mathcal{M}_p(\mathbb{R})\) by (Kraus and Panaretos 2012; Cardot and Godichon-Baggioni 2015):

\[ G_m(M) := \mathbb{E} \left[ \left\| (X - m)(X - m)^\top - M \right\|_{F} - \left\| (X - m)(X - m)^\top \right\|_{F} \right], \] where \(\| . \|_{F}\) denotes the Frobenius norm for matrices. The median covariance matrix can be computed using one of two available algorithms, implemented respectively in the functions `WeiszfeldMedianCovariance` and `ASGMedianCovariance`. When the distribution of \(X\) is symmetric, the median covariance matrix \(V\) and the usual covariance matrix \(\Sigma\) are not equal, but they share the same eigenvectors (see Proposition 2 in (Kraus and Panaretos 2012)). Moreover, denoting by \(\delta\) the vector containing the eigenvalues of \(V\) and \(\lambda\) the one for \(\Sigma\), one has (Proposition 2 in (Kraus and Panaretos 2012)):

\[ \delta = \frac{\mathbb{E} \left[ \frac{\lambda \odot U^{\otimes 2}}{h\left(\lambda,\delta,U \right)} \right] }{ \mathbb{E} \left[ \frac{1}{h\left(\lambda,\delta,U \right)} \right]}. \]

where \(U \sim \Sigma^{-1/2}(X-m)\), \(\lambda \odot U\) denotes the Hadamard product of the vectors \(\lambda\) and \(U\), and \(U^{\otimes 2} = U \odot U\). Additionally,

\[\begin{equation}\label{def::h} h(\lambda, \delta, U) := \sum_k (\delta[k] - \lambda[k] U[k]^2)^2 + \sum_{i \neq j} U[i]^2 U[j]^2 {\lambda[i]} {\lambda[j]}. \end{equation}\]

Then, one can use this relation to obtain robust estimates of the eigenvalues of the variance, and thereafter a robust estimate of the variance \(\Sigma\).

# 2 Robust estimates of the eigenvalues of the variance

The relation between \(\delta\) and \(\lambda\) given above allows us to reformulate the problem of finding \(\lambda\) as the search for the zero of a function. More precisely, denoting \(A = \operatorname{diag}(U[1]^2, \ldots, U[p]^2)\), the objective is to determine \(\lambda\) such that:

\[ \mathbb{E}\left[ \frac{A \lambda - \delta}{ {h(\lambda, \delta, U)}} \right] = 0. \]

Let us also suppose from now that one can generate i.i.d. copies \(U_{1}, \ldots , U_{N}\) of \(U\). For instance, in the case where \(X\) is a Gaussian vector, it means that \(U\) is a standardized Gaussian vector, and one can easily simulate \(N\) copies. Then, one can consider the empirical version of the above equation, and three algorithms for estimating the eigenvalues of \(\Sigma\) are available (see (Godichon-Baggioni and Robin 2022)): a Robbins–Monro algorithm, a gradient-type algorithm, and a fixed-point algorithm.

## 2.1 Robbins–Monro algorithm for estimating robustly the eigenvalues of the variance

Of course, one usually does not have access to the eigenvalues of the Median Covariation Matrix (MCM), nor to the MCM itself. Then, one can estimate the MCM with the help of Algorithm `weiszfeldMedianCovariance` or `ASGMedianCovariance`, before approximating its eigenvalues, to obtain an estimate \(\hat{\delta}\) of the vector of eigenvalues. The Robbins–Monro procedure (Robbins and Monro 1951; Godichon-Baggioni and Robin 2022) for estimating the eigenvalues of the variance can then be described step by step as follows:

1. Input:

   - \(U_1, \dots, U_N \in \mathbb{R}^p\),
   - initial value \(\lambda_{0} = \overline{\lambda}_{0} \in \mathbb{R}^p\),
   - Estimated eigenvalues of the MCM \(\hat{\delta} \in \mathbb{R}^p\),
   - Define constants:
     \(c_\gamma > 0\), \(\gamma \in (0.5, 1)\), \(n_0 \ge 0\), \(\omega \ge 0\)
   - Set \(S_0 = 1\)

2. Repeat for \(n = 0, \dots, N - 1\):

   1. Compute
      \(A_{n+1} = \mathrm{diag}(U_{n+1}[1]^2, \dots, U_{n+1}[p]^2)\)
   2. Compute the step size
      \(\gamma_{n+1} = c_\gamma (n + 1 + n_0)^{-\gamma}\)
   3. Update the current estimate
      \[ \lambda_{n+1} = \lambda_n - \gamma_{n+1} \frac{A_{n+1}\lambda_n - \hat{\delta}} { \sqrt{h(\lambda_n, \hat{\delta}, U_{n+1})}} \]
   4. Compute the weight
      \(w_{n+1} = \log(n + 1)^{\omega}\)
   5. Update the cumulative weight
      \(S_{n+1} = S_n + w_{n+1}\)
   6. Update the averaged estimate
      \[ \overline{\lambda}_{n+1} = \overline{\lambda}_n + \frac{w_{n+1}}{S_{n+1}} (\lambda_{n+1} - \overline{\lambda}_n) \]

3. Return the final averaged estimate \(\overline{\lambda}_N\)

The algorithm is implemented in the function `RobbinsMC`, which takes the following arguments:

- `U`: a matrix of size \(N \times p\) corresponding to a sample of \(\Sigma^{-1/2}(X - \mu)\). Each row represents one observation.

- `delta`: a vector of length `p` containing the eigenvalues of the median covariation matrix (MCM).

- `init`: initial value of the eigenvalues of the variance–covariance matrix (by default equal to `delta`).

- `gamma`: Robbins–Monro parameter (default: `0.75`).

- `c`: Robbins–Monro parameter (default: `2`).

- `w`: Robbins–Monro parameter (default: `2`).

- `epsilon`: stopping criterion. The algorithm stops when the difference between two iterations is less than `epsilon` (default: `1e-8`).

- `out_index`: indexes of iterations for which we want to store the estimates (default: number of observations `N`).

The function returns a list containing:
- `lambda`: the estimated eigenvalues,
- `iter`: the number of iterations performed,
- `vp`: the values of intermediate estimates (if requested).

### 2.1.1 Example of use

    # Define the true covariance matrix
    Sigma <- diag(2:1)
    N <- 5000

    # Generate Gaussian data
    X <- mvtnorm::rmvnorm(N, sigma = Sigma)

    # Compute the geometric median
    med <- STARRS::WeiszfeldMedian(X)

    # Compute the Median Covariation Matrix (MCM)
    MCM <- STARRS::WeiszfeldMedianCovariance(X, median_est = med)

    # Extract its eigenvalues
    delta <- eigen(MCM)$values
    p <- length(delta)

    # Generate U = Sigma^{-1/2}(X - mu)
    U <- matrix(rnorm(N * p), ncol = p)

    # Run the Robbins–Monro estimation
    lambda <- STARRS::robbinsMC(U = U, delta = delta, init = delta)

    # Display the estimated eigenvalues
    lambda
    #>           [,1]
    #> [1,] 1.8421897
    #> [2,] 0.9256411

## 2.2 Fixed-point algorithm for estimating the eigenvalues of the variance

Given an estimate \(\hat{\delta}\) of the vector of the eigenvalues of the Median Covariance Matrix, the fixed-point algorithm (Godichon-Baggioni and Robin 2022) provides a robust method to estimate the eigenvalues of the variance–covariance matrix.

The algorithm proceeds as follows:

1. Input:
   - Observations \(U_1, \dots, U_N \in \mathbb{R}^p\)
   - Initial value \(\lambda^{(0)} \in \mathbb{R}^p\)
   - Estimated eigenvalues of the MCM \(\hat{\delta} \in \mathbb{R}^p\)
   - Tolerance \(\varepsilon > 0\)
   - Maximum number of iterations \(T_{\max}\)
2. Iterate until convergence or maximum iterations:
   1. Initialize \(\text{num} = 0\) and \(\text{den} = 0\)
   2. For each \(i = 1, \dots, N\):
      - Compute
        \[ w_i =1/ \sqrt{h(\hat{\delta}, \lambda^{(k)}, U_i)}\]
      - Update accumulators:
        \(\text{num} = \text{num} + w_i\)
        \(\text{den} = \text{den} + U_i^{\otimes 2} * w_i\) (where \(\otimes\) is the elementwise product)
   3. Update
      \(\lambda^{(k+1)} = \hat{\delta} * (\text{num} \oslash \text{den})\)
      (where \(\oslash\) denotes coordinate-wise division)
   4. If \(\|\lambda^{(k+1)} - \lambda^{(k)}\| \leq \varepsilon\), stop.
3. Return: the final estimate \(\lambda^{(k+1)}\).

The algorithm is implemented in the function `FixMC`, which estimates the eigenvalues of the variance–covariance matrix based on a fixed-point procedure (Godichon-Baggioni and Robin 2022).

The function takes the following arguments:

- `U`: a matrix of size \(N \times p\) corresponding to \(\Sigma^{-1/2}(X - \mu)\). The rows are the observations.

- `delta`: a vector of size `p` containing the eigenvalues of the Median Covariation Matrix (MCM).

- `init`: the initial value of the vector of eigenvalues of the variance–covariance matrix (default: `delta`).

- `niter`: maximum number of iterations for the algorithm (default: `10`).

- `epsilon`: stopping criterion. The algorithm stops when the difference between two iterations is less than `epsilon` (default: `1e-8`).

- `out_index`: indexes of the iterations for which the estimates should be output (default: last iteration).

The function returns a list containing:
- `vp`: the estimated eigenvalues,
- `iter`: the number of iterations performed,
- `vplist`: the values of intermediate estimates (if requested).

### 2.2.1 Example of use

    # Define the true covariance matrix
    Sigma <- diag(2:1)
    N <- 5000

    # Generate Gaussian data
    X <- mvtnorm::rmvnorm(N, sigma = Sigma)

    # Compute the geometric median
    med <- STARRS::WeiszfeldMedian(X)

    # Compute the Median Covariation Matrix (MCM)
    MCM <- STARRS::WeiszfeldMedianCovariance(X, median_est = med)

    # Extract its eigenvalues
    delta <- eigen(MCM)$values
    p <- length(delta)

    # Generate standardized observations U = Sigma^{-1/2}(X - mu)
    U <- matrix(rnorm(N * p), ncol = p)

    # Run the fixed-point estimation
    lambda <- STARRS::fixMC(U = U, delta = delta, init = delta)

    # Display the estimated eigenvalues
    lambda
    #>           [,1]
    #> [1,] 1.8870428
    #> [2,] 0.8969555

## 2.3 Gradient-type algorithm

Given an estimate \(\hat{\delta}\) of the vector of the eigenvalues of the Median Covariance Matrix, the gradient-type algorithm (Godichon-Baggioni and Robin 2022) estimates the eigenvalues of the variance–covariance matrix as follows:

1. Input: sample \(U_1, \dots, U_N \in \mathbb{R}^p\),
   initial value \(\lambda^{(0)} \in \mathbb{R}^p\),
   Estimated eigenvalues of the MCM \(\hat{\delta} \in \mathbb{R}^p\),
   step-size sequence \(\{\eta_k\}_{k \geq 0}\),
   tolerance \(\varepsilon > 0\),
   and maximum number of iterations \(T_{\max}\).

2. For each iteration \(k = 0, 1, \dots, T_{\max}-1\):

   1. Initialize the gradient \(\leftarrow 0\) (vector of length \(p\)).
   2. For each observation \(i = 1, \dots, N\):
      • Compute \(w_i \leftarrow \sqrt{h(\hat{\delta}, \lambda^{(k)}, U_i)}\)
      • Compute \(g_i \leftarrow \left(U_i^{\otimes 2} \otimes \lambda^{(k)} - \hat{\delta}\right) / w_i\)
      • Update gradient \(\leftarrow\) gradient \(+ g_i\)
   3. Update \(\lambda^{(k+1)} \leftarrow \lambda^{(k)} - \eta_k \cdot \text{gradient}\)
   4. If \(\|\lambda^{(k+1)} - \lambda^{(k)}\| \leq \varepsilon\), stop.

3. Output: \(\lambda^{(k+1)}\)

The algorithm is implemented in the function `gradMC` with the following arguments:

- `U`: a matrix of size \(N \times p\) corresponding to \(\Sigma^{-1/2}(X - \mu)\). The rows are the observations.
- `delta`: vector of length `p` containing the eigenvalues of the MCM.
- `init`: initial vector of eigenvalues of the variance–covariance matrix (default: `delta`).
- `niter`: maximum number of iterations for the gradient descent algorithm (default: `10`).
- `epsilon`: stopping criterion — stops when the difference between two iterations is less than `1e-8`.
- `step`: step size of the gradient descent (default: `1`).
- `out_index`: indexes of iterations for which the estimates are returned (default: last iteration).

The function returns a list containing:
- `vp`: the estimated eigenvalues,
- `iter`: the number of iterations performed,
- `vplist`: the values of intermediate estimates (if requested).

### 2.3.1 Example of use

    # Load required packages
    library(mvtnorm)
    #>
    #> Attaching package: 'mvtnorm'
    #> The following object is masked from 'package:mclust':
    #>
    #>     dmvnorm
    library(STARRS)

    # Define the true covariance matrix
    Sigma <- diag(2:1)
    N <- 5000

    # Generate Gaussian data
    X <- mvtnorm::rmvnorm(N, sigma = Sigma)

    # Compute the geometric median
    med <- STARRS::WeiszfeldMedian(X)

    # Compute the Median Covariation Matrix (MCM)
    MCM <- STARRS::WeiszfeldMedianCovariance(X, median_est = med)

    # Extract eigenvalues of the MCM
    delta <- eigen(MCM)$values
    p <- length(delta)

    # Generate standardized observations
    U <- matrix(rnorm(N * p), ncol = p)

    # Run the gradient algorithm
    lambda <- STARRS::gradMC(U = U, delta = delta, init = delta)

    # Display the estimated eigenvalues
    lambda
    #>           [,1]
    #> [1,] 1.7178355
    #> [2,] 0.8730577

## 2.4 Comparison of the methods

In order to compare the different methods, we consider a random variable \(X \sim \mathcal{N}\left( 0 , \Sigma \right)\) where \(\Sigma = \text{Diag} \left( 1 , \sqrt{2} , \ldots , \sqrt{d} \right)\) and \(d = 20\). To obtain an accurate estimate of the Median Covariance Matrix, we first generate \(n = 10^{6}\) independent realizations of \(X\) and apply the Weiszfeld algorithm to estimate the MCM. The eigenvalues of this estimated MCM are computed using the R function `eigen`. Then, to assess and compare the three algorithms, we generate \(100\) replications of \(U \sim \mathcal{N}(0, I_d)\) for several sample sizes \(N \in \{100, 200, 2000, 5000, 10000\}\).

For each method, we compute:
1. The squared error between the estimated and true eigenvalues of the covariance. 2. The computational time (in seconds).

Figure below shows the comparison of squared errors and computation times. We observe that the Robbins–Monro method is slightly less accurate than the others, but its computation time is divided by about 15. Hence, it appears to be the preferred choice in large-scale or streaming contexts.

Figure 2.1: Comparison of the squared errors and computational time between the three methods (Robbins–Monro, Gradient, and Fixed-point) as a function of the sample size.

# 3 Robust estimation of the variance

## 3.1 Offline Robust Variance Estimation

In all the sequel, we consider that \(X\) is a random vector of mean \(\mu\) and variance \(\Sigma\), and we assume that it is possible to simulate according to \(U = \Sigma^{-1/2} (X - \mu)\). Given a sample \(X_1, \ldots, X_N\), the robust estimate of the variance consists in obtaining an estimate of the median and the MCM with the help of the Weiszfeld algorithms given by Algorithms `WeiszfeldMedian` and `WeiszfeldMedianCovariance`, before approximating the eigenvalues (resp. eigenvectors) of the MCM with the help of the R function `eigen` and we denote by \(\hat{\delta}\) (resp. \(\hat{P}\)) these approximations. Then one can generate \(n\) independent realizations \(U_1, \ldots, U_n\) with the same law as \(U\). Then, one can choose one of the three methods (Robbins-Monro given by Algorithm `robbinsMC`, gradient given by Algorithm `gradMC` or fix point given by Algorithm `fixMC`) for obtaining an estimate \(\hat{\lambda}\) of the eigenvalues of the covariance. Then, one can obtain a robust estimate of the variance with

\[ \hat{\Sigma} = \hat{P} \, \text{Diag}(\hat{\lambda}) \, \hat{P}^{T}. \]

This is given by the following Algorithm:

1. Input: sample \(X_1, \dots, X_N\) Number of Monte Carlo iterations \(n_{MC}\)

2. Compute the robust median:

   - `mu_hat <- WeiszfeldMedian(X)`

3. Compute the Median Covariation Matrix (MCM):

   - `Gamma_hat <- WeiszfeldMedianCovariance(X, median_est = mu_hat)`

4. Eigen-decomposition of the MCM:

   - `(delta_hat, P_hat) <- eigen(Gamma_hat)`
     - `delta_hat`: eigenvalues
     - `P_hat`: eigenvectors

5. Generate Monte Carlo samples \(\mathbf{U} = \left( U_1, \dots, U_{n_{MC}} \right)\) with \(U_{i} \sim U\).

6. Estimate eigenvalues of the variance using one of the three methods:

   - `lambda_hat <- RobbinsMC(\mathbf{U}, delta_hat)` or
   - `lambda_hat <- gradMC(\mathbf{U}, delta_hat)` or
   - `lambda_hat <- fixMC(\mathbf{U}, delta_hat)`

7. Reconstruct the robust variance estimate:
   \[ \hat{\Sigma} = P_{\text{hat}} \cdot \text{Diag}(\lambda_{\text{hat}}) \cdot P_{\text{hat}}^T \]

8. Output:

   - `Sigma_hat` : estimated covariance
   - `mu_hat` : median
   - `Gamma_hat` : MCM
   - `delta_hat`, `P_hat` : eigenvalues and eigenvectors of MCM

### 3.1.1 Example of use

    library(mvtnorm)
    library(STARRS)

    # Generate sample data
    N <- 5000
    p <- 3
    Sigma <- diag(1:p)
    X <- mvtnorm::rmvnorm(N, sigma = Sigma)

    # Compute offline robust variance
    Variance <- STARRS::onlineRobustVariance(X)

    # Display the estimated covariance
    Variance$variance
    #>             [,1]        [,2]        [,3]
    #> [1,]  0.84967546 -0.02096993 -0.06859078
    #> [2,] -0.02096993  1.80644317  0.01898015
    #> [3,] -0.06859078  0.01898015  2.58218512

## 3.2 Online robust estimation of the variance

The online robust estimation of the variance, or streaming approach, consists in simultaneously estimating both the median and the Median Covariance Matrix (MCM) using a streaming gradient algorithm (Godichon-Baggioni et al. 2023). More precisely, at each gradient step (used to estimate the median and the MCM), the algorithm approximates the eigenvalues and eigenvectors of the MCM. Starting from these estimates of the MCM’s eigenvalues and from the current estimate of the covariance’s eigenvalues, several iterations of the Robbins–Monro algorithm (Algorithm `robbinsMC``) are then performed to obtain a new estimate of the covariance eigenvalues. The covariance matrix can then be reconstructed as in the offline case.

This is given by the following Algorithm:

1. Input:
   - Streaming data blocks \(\left\{ \{ X_{n,j} \}_{j=1}^{k} \right\}_{n=1}^{N} \subset \mathbb{R}^p\)
   - Step sizes \(\gamma_n\)
   - Number of Monte Carlo iterations \(n_{MC}\)
   - Initial median \(m_0 = \overline{m}_{0}\)
   - Initial MCM \(V_0 = \bar{V}_0\)
   - Initial eigenvalues \(\lambda_0\)
2. Iterate for each block \(n = 0, \dots, N-1\):
   1. Retrieve next block of \(k\) observations: `X_block = {X_{n+1,j}}`
   2. Update median: \[ m_{n+1} = m_n + \gamma_{n+1} \frac{1}{k} \sum_{j=1}^{k} \frac{X_{n+1,j} - m_n}{\| X_{n+1,j} - m_n \|}, \] \[ \bar{m}_{n+1} = \bar{m}_n + \frac{1}{n+2} ( m_{n+1} - \bar{m}_n ) \]
   3. Update MCM (matrix of second moments): \[ V_{n+1} = V_n + \gamma_{n+1} \frac{1}{k} \sum_{j=1}^{k} \frac{(X_{n+1,j} - \bar{m}_n)(X_{n+1,j} - \bar{m}_n)^\top - V_n}{\| \cdot \|_F}, \] \[ \bar{V}_{n+1} = \bar{V}_n + \frac{1}{n+2} ( V_{n+1} - \bar{V}_n ) \]
   4. Eigen-decomposition of averaged MCM: \[ \delta_{n+1} = \text{eigenvalues}(\bar{V}_{n+1}), \quad P_{n+1} = \text{eigenvectors}(\bar{V}_{n+1}) \]
   5. Monte Carlo update of covariance eigenvalues:
      - Generate \(\mathbf{U} = \left( U_1, \dots, U_{n_{MC}} \right)\) where \(U_{i} \sim U\)
      - Update \[ \lambda_{n+1} = \text{robbinsMC}(\mathbf{U}, \delta_{n+1}, \lambda_n) \]
   6. Reconstruct covariance: \[ \Sigma_{n+1} = P_{n+1} \, \text{diag}(\lambda_{n+1}) \, P_{n+1}^\top \]
3. Return:
   - \(\Sigma_N\): estimated covariance
   - \(\bar{m}_N\): estimated median
   - \(\bar{V}_N\): estimated MCM
   - \(\lambda_N\): estimated eigenvalues

### 3.2.1 Example of use

    library(mvtnorm)
    library(STARRS)

    # Generate sample data
    N <- 5000
    p <- 3
    Sigma <- diag(1:p)
    X <- mvtnorm::rmvnorm(N, sigma = Sigma)

    # Compute online robust variance
    Variance <- STARRS::onlineRobustVariance(X)

    # Display the estimated covariance
    Variance$variance
    #>             [,1]        [,2]        [,3]
    #> [1,]  0.88779232 -0.02594460  0.01254807
    #> [2,] -0.02594460  1.73499517 -0.04026758
    #> [3,]  0.01254807 -0.04026758  2.80287366

## 3.3 Comparison of the methods

In order to compare the different methods, we consider a random variable \(X \sim \mathcal{N}(0, \Sigma)\) where \(\Sigma = \text{Diag}(1, \sqrt{2}, \ldots, \sqrt{d})\) and \(d = 20\). We generate 100 samples (of \(U \sim \mathcal{N}(0, I_d)\)) of sizes 1000, 2000, 5000, 10000 before considering the squared error of the estimation of the covariance as well as the computational time for each method.

Figure 3.1: Comparison of the different methods for estimating robustly the variance

# 4 References

Cardot, Hervé, and Antoine Godichon-Baggioni. 2015. “Fast Estimation of the Median Covariation Matrix with Application to Online Robust Principal Components Analysis.” TEST, 1–20.
Godichon-Baggioni, Antoine, Nicklas Werge, and Olivier Wintenberger. 2023. “Non-Asymptotic Analysis of Stochastic Approximation Algorithms for Streaming Data.” ESAIM: Probability and Statistics 27: 482–514.
Godichon-Baggioni, A., and S. Robin. 2022. “A Robust Model-Based Clustering Based on the Geometric Median and the Median Covariation Matrix.” Statistics and Computing 34 (1): 55.
Haldane, J. B. S. 1948. “Note on the median of a multivariate distribution.” Biometrika 35 (3-4): 414–17. https://doi.org/10.1093/biomet/35.3-4.414.
Kraus, David, and Victor M. Panaretos. 2012. “Dispersion Operators and Resistant Second-Order Functional Data Analysis.” Biometrika 99: 813–32.
Robbins, Herbert, and Sutton Monro. 1951. “A Stochastic Approximation Method.” The Annals of Mathematical Statistics, 400–407.