LeJEPA: Provable and Scalable Self-Supervised Learning Without the Heuristics — part 9/10

Cleaned-text code points [117672, 132698). Bibliography omitted here; source.md includes it.

A direct calculation shows Fix t ∈ ℝ d t\in\mathbb{R}^{d} and abbreviate Z j ≔ e i ​ t ⊤ ​ X j Z_{j}\coloneqq e^{\mathrm{i}t^{\top}X_{j}} , so that ϕ n ​ ( t ) = 1 n ​ ∑ j = 1 n Z j \phi_{n}(t)=\frac{1}{n}\sum_{j=1}^{n}Z_{j} . Note that | Z j | = 1 |Z_{j}|=1 almost surely (since t ⊤ ​ X j ∈ ℝ t^{\top}X_{j}\in\mathbb{R} ), and 𝔼 ⁡ [ Z j ] = ϕ θ ​ ( t ) \mathbb{E}[Z_{j}]=\phi_{\theta}(t) for all j j . We start from the algebraic identity | ϕ n ​ ( t ) − ψ ⁡ ( t ) | 2 = ϕ n ​ ( t ) ​ ϕ n ​ ( t ) ¯ − ψ ⁡ ( t ) ​ ϕ n ​ ( t ) ¯ − ψ ⁡ ( t ) ¯ ​ ϕ n ​ ( t ) + | ψ ⁡ ( t ) | 2 . \big|\phi_{n}(t)-\psi(t)\big|^{2}=\phi_{n}(t)\overline{\phi_{n}(t)}-\psi(t)\overline{\phi_{n}(t)}-\overline{\psi(t)}\phi_{n}(t)+\big|\psi(t)\big|^{2}. Taking expectations term by term gives 𝔼 ⁡ [ | ϕ n − ψ | 2 ] = \displaystyle\mathbb{E}\left[\big|\phi_{n}-\psi\big|^{2}\right]= 𝔼 ⁡ [ | ϕ n | 2 ] − ψ ​ 𝔼 ​ [ ϕ n ¯ ] − ψ ¯ ​ 𝔼 ​ [ ϕ n ] + | ψ | 2 , \displaystyle\mathbb{E}\left[|\phi_{n}|^{2}\right]-\psi\mathbb{E}\left[\overline{\phi_{n}}\right]-\overline{\psi}\mathbb{E}\left[\phi_{n}\right]+|\psi|^{2}, (28) = \displaystyle= 𝔼 ⁡ [ | ϕ n | 2 ] − ψ ​ 𝔼 ⁡ [ ϕ n ] ¯ − ψ ¯ ​ 1 n ​ ∑ j = 1 n 𝔼 ⁡ [ Z j ] + | ψ | 2 , \displaystyle\mathbb{E}\left[|\phi_{n}|^{2}\right]-\psi\overline{\mathbb{E}[\phi_{n}]}-\overline{\psi}\frac{1}{n}\sum_{j=1}^{n}\mathbb{E}[Z_{j}]+|\psi|^{2}, (29) = \displaystyle= 𝔼 ⁡ [ | ϕ n | 2 ] − ψ ​ ϕ θ ¯ − ψ ¯ ​ ϕ θ + | ψ | 2 , \displaystyle\mathbb{E}\left[|\phi_{n}|^{2}\right]-\psi\overline{\phi_{\theta}}-\overline{\psi}\phi_{\theta}+|\psi|^{2}, (30) = \displaystyle= 𝔼 ⁡ [ | ϕ n | 2 ] − 2 ​ R ​ e ​ ( ψ ¯ ​ ϕ θ ) + | ψ | 2 , \displaystyle\mathbb{E}\left[|\phi_{n}|^{2}\right]-2\mathrm{Re}\big(\overline{\psi}\phi_{\theta}\big)+|\psi|^{2}, (31) = \displaystyle= 𝔼 ⁡ [ | 1 n ​ ∑ j = 1 n Z j | 2 ] − 2 ​ R ​ e ​ ( ψ ¯ ​ ϕ θ ) + | ψ | 2 , \displaystyle\mathbb{E}\left[\left|\frac{1}{n}\sum_{j=1}^{n}Z_{j}\right|^{2}\right]-2\mathrm{Re}\big(\overline{\psi}\phi_{\theta}\big)+|\psi|^{2}, (32) = \displaystyle= 1 n 2 ​ ∑ j = 1 n ∑ l = 1 n 𝔼 ⁡ [ Z j ​ Z l ¯ ] − 2 ​ R ​ e ​ ( ψ ¯ ​ ϕ θ ) + | ψ | 2 , \displaystyle\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{l=1}^{n}\mathbb{E}\left[Z_{j}\overline{Z_{l}}\right]-2\mathrm{Re}\big(\overline{\psi}\phi_{\theta}\big)+|\psi|^{2}, (33) Since the Z j Z_{j} are i.i.d., 𝔼 ⁡ [ Z j ​ Z l ¯ ] = { 𝔼 ⁡ [ | Z 1 | 2 ] = 1 , if ​ j = l , 𝔼 ⁡ [ Z j ] ​ 𝔼 ⁡ [ Z l ] ¯ = ϕ θ ​ ϕ θ ¯ = | ϕ θ | 2 , if ​ j ≠ l , \mathbb{E}\left[Z_{j}\overline{Z_{l}}\right]=\begin{cases}\mathbb{E}\left[|Z_{1}|^{2}\right]=1,&\text{if }j=l,\\[4.0pt] \mathbb{E}[Z_{j}]\overline{\mathbb{E}[Z_{l}]}=\phi_{\theta}\overline{\phi_{\theta}}=|\phi_{\theta}|^{2},&\text{if }j\neq l,\end{cases} hence 𝔼 ⁡ [ | ϕ n | 2 ] = \displaystyle\mathbb{E}\left[|\phi_{n}|^{2}\right]= 1 n 2 ​ ( n + n ⁡ ( n − 1 ) ​ | ϕ θ | 2 ) \displaystyle\frac{1}{n^{2}}\Big(n+n(n-1)|\phi_{\theta}|^{2}\Big) = \displaystyle= 1 n + ( 1 − 1 n ) ​ | ϕ θ | 2 \displaystyle\frac{1}{n}+\left(1-\frac{1}{n}\right)|\phi_{\theta}|^{2} = \displaystyle= | ϕ θ | 2 + 1 − | ϕ θ | 2 n \displaystyle|\phi_{\theta}|^{2}+\frac{1-|\phi_{\theta}|^{2}}{n} Plugging these, we obtain 𝔼 ⁡ [ | ϕ n − ψ | 2 ] \displaystyle\mathbb{E}\left[\big|\phi_{n}-\psi\big|^{2}\right] = ( | ϕ θ | 2 + 1 − | ϕ θ | 2 n ) − 2 ​ R ​ e ​ ( ψ ¯ ​ ϕ θ ) + | ψ | 2 \displaystyle=\left(|\phi_{\theta}|^{2}+\frac{1-|\phi_{\theta}|^{2}}{n}\right)-2\mathrm{Re}\big(\overline{\psi}\phi_{\theta}\big)+|\psi|^{2} = ( | ϕ θ | 2 − 2 ​ R ​ e ​ ( ψ ¯ ​ ϕ θ ) + | ψ | 2 ) + 1 − | ϕ θ | 2 n \displaystyle=\big(|\phi_{\theta}|^{2}-2\mathrm{Re}\big(\overline{\psi}\phi_{\theta}\big)+|\psi|^{2}\big)+\frac{1-|\phi_{\theta}|^{2}}{n} = | ϕ θ − ψ | 2 + 1 − | ϕ θ | 2 n . \displaystyle=\big|\phi_{\theta}-\psi\big|^{2}+\frac{1-|\phi_{\theta}|^{2}}{n}. Under Dominated convergence, 𝔼 ⁡ [ ∇ θ D n ​ ( t ) ] = ∇ θ 𝔼 ​ [ D n ​ ( t ) ] \mathbb{E}[\nabla_{\theta}D_{n}(t)]=\nabla_{\theta}\mathbb{E}[D_{n}(t)] , hence 𝔼 ⁡ [ ∇ θ D n ​ ( t ) ] = ∇ θ | ϕ θ ​ ( t ) − ψ ⁡ ( t ) | 2 + ∇ θ 1 − | ϕ θ ​ ( t ) | 2 n , \mathbb{E}\left[\nabla_{\theta}D_{n}(t)\right]=\nabla_{\theta}\big|\phi_{\theta}(t)-\psi(t)\big|^{2}+\nabla_{\theta}\frac{1-|\phi_{\theta}(t)|^{2}}{n}, concluding the proof.

In practice one replaces ∫ ℝ w ​ ( t ) ​ ( ⋅ ) ​ 𝑑 t \int_{\mathbb{R}}w(t)(\cdot)dt by a deterministic quadrature on a uniform grid t k ∈ [ − T , T ] t_{k}\in[-T,T] with weights ω k \omega_{k} (e.g. trapezoidal rule) and a Gaussian window w ⁡ ( t ) = e − α ​ t 2 w(t)=e^{-\alpha t^{2}} . All statements above remain valid with the integral replaced by ∑ k ω k ​ ( ⋅ ) \sum_{k}\omega_{k}(\cdot) : L ⁡ ( θ ) ≈ ∑ k ω k ​ | ϕ θ ​ ( t k ) − ψ ⁡ ( t k ) | 2 , L ^ n ​ ( θ ) ≈ ∑ k ω k ​ | ϕ n ​ ( t k ) − ψ ⁡ ( t k ) | 2 , L(\theta)\approx\sum_{k}\omega_{k}\big|\phi_{\theta}(t_{k})-\psi(t_{k})\big|^{2},\quad\widehat{L}_{n}(\theta)\approx\sum_{k}\omega_{k}\big|\phi_{n}(t_{k})-\psi(t_{k})\big|^{2}, and the bias term becomes Bias ( θ ) = − 1 n ∑ k ω k ∇ θ | ϕ θ ( t k ) | 2 . \text{Bias}(\theta)=-\frac{1}{n}\sum_{k}\omega_{k}\nabla_{\theta}\big|\phi_{\theta}(t_{k})\big|^{2}. Since the grid and weights are deterministic, they do not affect unbiasedness with respect to sampling; they only introduce a deterministic approximation error to the target functional L ⁡ ( θ ) L(\theta) .

∎

### B.14 Proof of VICReg’s Recovery

###### Proof.

We prove this result in two parts.

##### Part I: 𝔼 ⁡ [ 𝐗 ] = 𝟎 \mathbb{E}[\mathbf{X}]=\mathbf{0}

Given that 𝔼 ⁡ [ ⟨ 𝐗 , 𝐚 ⟩ ] = 0 \mathbb{E}[\langle\mathbf{X},\mathbf{a}\rangle]=0 for all unit vectors 𝐚 \mathbf{a} , and noting that ⟨ 𝐗 , 𝐚 ⟩ = 𝐚 T ​ 𝐗 \langle\mathbf{X},\mathbf{a}\rangle=\mathbf{a}^{T}\mathbf{X} , we have: 𝔼 ⁡ [ 𝐚 T ​ 𝐗 ] = 0 for all ​ 𝐚 ∈ ℝ d ​ with ​ ‖ 𝐚 ‖ = 1 \mathbb{E}[\mathbf{a}^{T}\mathbf{X}]=0\quad\text{for all }\mathbf{a}\in\mathbb{R}^{d}\text{ with }\|\mathbf{a}\|=1 (35) By linearity of expectation: 𝐚 T ​ 𝔼 ​ [ 𝐗 ] = 0 for all unit vectors ​ 𝐚 \mathbf{a}^{T}\mathbb{E}[\mathbf{X}]=0\quad\text{for all unit vectors }\mathbf{a} (36) Let 𝝁 = 𝔼 ⁡ [ 𝐗 ] \bm{\mu}=\mathbb{E}[\mathbf{X}] . We claim that 𝝁 = 𝟎 \bm{\mu}=\mathbf{0} . Suppose, for the sake of contradiction, that 𝝁 ≠ 𝟎 \bm{\mu}\neq\mathbf{0} . Then ‖ 𝝁 ‖ 2 > 0 \|\bm{\mu}\|_{2}>0 . Define the unit vector: 𝐚 ∗ = 𝝁 ‖ 𝝁 ‖ 2 \mathbf{a}^{*}=\frac{\bm{\mu}}{\|\bm{\mu}\|_{2}} (37) Since 𝐚 ∗ \mathbf{a}^{*} is a unit vector, equation ( 35 ) implies: ( 𝐚 ∗ ) T ​ 𝝁 = 0 (\mathbf{a}^{*})^{T}\bm{\mu}=0 (38) However, substituting the definition of 𝐚 ∗ \mathbf{a}^{*} : ( 𝐚 ∗ ) T ​ 𝝁 = ( 𝝁 ‖ 𝝁 ‖ 2 ) T ​ 𝝁 = 𝝁 T ​ 𝝁 ‖ 𝝁 ‖ 2 = ‖ 𝝁 ‖ 2 2 ‖ 𝝁 ‖ 2 = ‖ 𝝁 ‖ 2 > 0 (\mathbf{a}^{*})^{T}\bm{\mu}=\left(\frac{\bm{\mu}}{\|\bm{\mu}\|_{2}}\right)^{T}\bm{\mu}=\frac{\bm{\mu}^{T}\bm{\mu}}{\|\bm{\mu}\|_{2}}=\frac{\|\bm{\mu}\|_{2}^{2}}{\|\bm{\mu}\|_{2}}=\|\bm{\mu}\|_{2}>0 (39) This contradiction establishes that 𝝁 = 𝟎 \bm{\mu}=\mathbf{0} .

##### Part II: Cov ⁡ ( 𝐗 ) = 𝐈 d \mathrm{Cov}(\mathbf{X})=\mathbf{I}_{d}

Since 𝔼 ⁡ [ 𝐗 ] = 𝟎 \mathbb{E}[\mathbf{X}]=\mathbf{0} , we have: Var ⁡ ( ⟨ 𝐗 , 𝐚 ⟩ ) = 𝔼 ⁡ [ ( ⟨ 𝐗 , 𝐚 ⟩ ) 2 ] = 𝔼 ⁡ [ ( 𝐚 T ​ 𝐗 ) 2 ] \mathrm{Var}(\langle\mathbf{X},\mathbf{a}\rangle)=\mathbb{E}[(\langle\mathbf{X},\mathbf{a}\rangle)^{2}]=\mathbb{E}[(\mathbf{a}^{T}\mathbf{X})^{2}] (40) Expanding the quadratic form: 𝔼 ⁡ [ ( 𝐚 T ​ 𝐗 ) 2 ] = 𝔼 ⁡ [ 𝐚 T ​ 𝐗𝐗 T ​ 𝐚 ] = 𝐚 T ​ 𝔼 ​ [ 𝐗𝐗 T ] ​ 𝐚 \mathbb{E}[(\mathbf{a}^{T}\mathbf{X})^{2}]=\mathbb{E}[\mathbf{a}^{T}\mathbf{X}\mathbf{X}^{T}\mathbf{a}]=\mathbf{a}^{T}\mathbb{E}[\mathbf{X}\mathbf{X}^{T}]\mathbf{a} (41) Since 𝔼 ⁡ [ 𝐗 ] = 𝟎 \mathbb{E}[\mathbf{X}]=\mathbf{0} , the covariance matrix is Cov ⁡ ( 𝐗 ) = 𝔼 ⁡ [ 𝐗𝐗 T ] \mathrm{Cov}(\mathbf{X})=\mathbb{E}[\mathbf{X}\mathbf{X}^{T}] . Let 𝚺 = Cov ⁡ ( 𝐗 ) \bm{\Sigma}=\mathrm{Cov}(\mathbf{X}) . The variance condition gives us: 𝐚 T ​ 𝚺 ​ 𝐚 = 1 for all unit vectors ​ 𝐚 \mathbf{a}^{T}\bm{\Sigma}\mathbf{a}=1\quad\text{for all unit vectors }\mathbf{a} (42) We now show that 𝚺 = 𝐈 d \bm{\Sigma}=\mathbf{I}_{d} . Step 1: Diagonal entries. For i ∈ { 1 , 2 , … , d } i\in\{1,2,\ldots,d\} , let 𝐞 i \mathbf{e}_{i} denote the i i -th standard basis vector. Setting 𝐚 = 𝐞 i \mathbf{a}=\mathbf{e}_{i} in equation ( 42 ): 𝐞 i T ​ 𝚺 ​ 𝐞 i = Σ i ​ i = 1 \mathbf{e}_{i}^{T}\bm{\Sigma}\mathbf{e}_{i}=\Sigma_{ii}=1 (43) Therefore, all diagonal entries of 𝚺 \bm{\Sigma} equal 1. Step 2: Off-diagonal entries. For distinct indices i , j ∈ { 1 , 2 , … , d } i,j\in\{1,2,\ldots,d\} , consider the unit vector: 𝐚 = 𝐞 i + 𝐞 j ‖ 𝐞 i + 𝐞 j ‖ 2 = 𝐞 i + 𝐞 j 2 \mathbf{a}=\frac{\mathbf{e}_{i}+\mathbf{e}_{j}}{\|\mathbf{e}_{i}+\mathbf{e}_{j}\|_{2}}=\frac{\mathbf{e}_{i}+\mathbf{e}_{j}}{\sqrt{2}} (44) Applying equation ( 42 ): 𝐚 T ​ 𝚺 ​ 𝐚 = 1 2 ​ ( 𝐞 i + 𝐞 j ) T ​ 𝚺 ​ ( 𝐞 i + 𝐞 j ) = 1 \mathbf{a}^{T}\bm{\Sigma}\mathbf{a}=\frac{1}{2}(\mathbf{e}_{i}+\mathbf{e}_{j})^{T}\bm{\Sigma}(\mathbf{e}_{i}+\mathbf{e}_{j})=1 (45) Expanding the quadratic form and using the symmetry of 𝚺 \bm{\Sigma} : 1 2 ​ ( 𝐞 i T ​ 𝚺 ​ 𝐞 i + 2 ​ 𝐞 i T ​ 𝚺 ​ 𝐞 j + 𝐞 j T ​ 𝚺 ​ 𝐞 j ) \displaystyle\frac{1}{2}(\mathbf{e}_{i}^{T}\bm{\Sigma}\mathbf{e}_{i}+2\mathbf{e}_{i}^{T}\bm{\Sigma}\mathbf{e}_{j}+\mathbf{e}_{j}^{T}\bm{\Sigma}\mathbf{e}_{j}) = 1 \displaystyle=1 (46) 1 2 ​ ( Σ i ​ i + 2 ​ Σ i ​ j + Σ j ​ j ) \displaystyle\frac{1}{2}(\Sigma_{ii}+2\Sigma_{ij}+\Sigma_{jj}) = 1 \displaystyle=1 (47) 1 2 ​ ( 1 + 2 ​ Σ i ​ j + 1 ) \displaystyle\frac{1}{2}(1+2\Sigma_{ij}+1) = 1 \displaystyle=1 (48) 1 + Σ i ​ j \displaystyle 1+\Sigma_{ij} = 1 \displaystyle=1 (49) Σ i ​ j \displaystyle\Sigma_{ij} = 0 \displaystyle=0 (50) Therefore, all off-diagonal entries of 𝚺 \bm{\Sigma} equal zero, establishing that 𝚺 = 𝐈 d \bm{\Sigma}=\mathbf{I}_{d} . ∎

## Appendix C Background

Foundation: The Linear Regression Model We start with the standard linear regression model: 𝐲 = 𝐗 ​ 𝜷 + 𝜺 \mathbf{y}=\mathbf{X}\bm{\beta}+\bm{\varepsilon} where: • 𝐲 = [ y 1 , y 2 , … , y n ] T ∈ ℝ n \mathbf{y}=[y_{1},y_{2},\ldots,y_{n}]^{T}\in\mathbb{R}^{n} is the response vector

• 𝐗 ∈ ℝ n × p \mathbf{X}\in\mathbb{R}^{n\times p} is the design matrix with 𝐗 i ​ j = x i ​ j \mathbf{X}_{ij}=x_{ij}

• 𝜷 = [ β 1 , β 2 , … , β p ] T ∈ ℝ p \bm{\beta}=[\beta_{1},\beta_{2},\ldots,\beta_{p}]^{T}\in\mathbb{R}^{p} is the parameter vector

• 𝜺 = [ ε 1 , ε 2 , … , ε n ] T ∼ 𝒩 ⁡ ( 𝟎 , σ 2 ​ 𝐈 n ) \bm{\varepsilon}=[\varepsilon_{1},\varepsilon_{2},\ldots,\varepsilon_{n}]^{T}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{n}) is the error vector

The error assumption means:

Step 1: Deriving the OLS Estimator To find the OLS estimator, we minimize the sum of squared residuals:

Expanding this quadratic form:

Taking the derivative with respect to 𝜷 \bm{\beta} :

Setting equal to zero and solving:

Assuming 𝐗 T ​ 𝐗 \mathbf{X}^{T}\mathbf{X} is invertible:

dimension=128, slices=10 dimension=128, slices=100 dimension=1024, slices=100

## Appendix D Details on Low-Discrepancy Sequences

Quasi-Monte Carlo (QMC) methods, such as the Sobol sequence, are widely used to generate low-discrepancy samples in the unit hypercube, providing improved uniformity over purely random sampling. To obtain samples uniformly distributed on the hypersphere, each QMC point is mapped to a standard normal vector via the inverse cumulative distribution function (CDF), and then projected onto the sphere by normalization. This approach leverages the rotational invariance of the multivariate normal distribution, ensuring that the resulting directions are uniformly distributed on the sphere’s surface. While the low-discrepancy property is not strictly preserved under this nonlinear mapping, the resulting samples are empirically more uniform than random samples and are standard in high-dimensional applications Marsaglia [1972] , Dick and Pillichshammer [2010] , Caflisch [1998] . 1: Number of points N N , dimension d d

## Appendix E Shapiro-Wilk Test

Let X1 < X2 < . . . < Xn denote an ordered random sample of size n from a standard normal distribution. Also, let m 5 (m1,m2,…,mn) be the vector of expected values of standard normal order statistics, and let V 5 (vij ) be the corresponding n 3 n covariance matrix, so that E ⁡ ( X i ) = m i and cov ⁡ ( X i , X j ) = v i ​ j , i , j = 1 , 2 , … , n E\left(X_{i}\right)=m_{i}\quad\text{ and }\quad\operatorname{cov}\left(X_{i},X_{j}\right)=v_{ij},\quad i,j=1,2,\ldots,n (52) The W test statistic Shapiro and Wilk [1965] for normality is then denoted by W = ( ∑ i = 1 n a i ​ Y i ) ∑ i = 1 n ( Y i − Y ¯ ) 2 = ( 𝐚𝐘 ) S 2 𝐚 ′ = ( a 1 , a 2 , … , a n ) = 𝐦𝐕 − 1 ( 𝐦𝐕 − 1 𝐕 − 1 𝐦 ) − 1 / 2 S 2 = ∑ i = 1 n ( Y i − Y ¯ ) 2 \begin{array}[]{l}W=\frac{\left(\sum_{i=1}^{n}a_{i}Y_{i}\right)}{\sum_{i=1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}}=\frac{(\mathbf{a}\mathbf{Y})}{S^{2}}\\ \mathbf{a}^{\prime}=\left(a_{1},a_{2},\ldots,a_{n}\right)=\mathbf{m}\mathbf{V}^{-1}\left(\mathbf{m}\mathbf{V}^{-1}\mathbf{V}^{-1}\mathbf{m}\right)^{-1/2}\\ \mathrm{S}^{2}=\sum_{i=1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}\end{array} (53)

Shapiro and Francia [1972] suggested replacing the covariance matrix V by the identity matrix I, because for large samples, the observations Yi may be treated as if they are independent (see Gupta [1952] ). Another asymptotic extension was suggested by Weisburg and Binham [1975] E ⁡ ( X i ) = m i ≈ Φ − 1 ​ ( i − 3 8 n + 1 4 ) i = 1 , 2 , … , n E\left(X_{i}\right)=m_{i}\approx\Phi^{-1}\left(\frac{i-\frac{3}{8}}{n+\frac{1}{4}}\right)\quad i=1,2,\ldots,n (54) building atop Elfving [1947] ’s approximation but using 3 / 8 3/8 instead of π / 8 \pi/8 .

Rahman and Govindarajulu [1997] proposed another variation using the approximation for the expected values of order statistics given by Blom [1958] and the approximations for the elements of the variance± covariance matrix given by Blom [1958] , Mosteller [2006] . These approximations are E ( X i ) = m i ≈ Φ − 1 ( i N + 1 ) , i = 1 , 2 , … , n E\left(X_{i}\right)=m_{i}\approx\Phi^{-1}\left(\frac{i}{N+1}\right),\quad i=1,2,\ldots,n (55) cov ⁡ ( X i , X j ) = v i ​ j ≈ p i ​ p j ( n + 2 ) ​ f ​ ( m i ) ​ f ​ ( m j ) , i , j = 1 , 2 , … , n \operatorname{cov}\left(X_{i},X_{j}\right)=v_{ij}\approx\frac{p_{i}p_{j}}{(n+2)f\left(m_{i}\right)f\left(m_{j}\right)},\quad i,j=1,2,\ldots,n (56) p i = i n + 1 p_{i}=\frac{i}{n+1} (57)

We know (see Hammersley and Morton [1954] , Plackett [1958] ) 𝐕 − 1 = ( n + 1 ) ​ ( n + 2 ) × ( 2 ​ ϕ 2 ​ ( m 1 ) − ϕ ⁡ ( m 1 ) ​ ϕ ​ ( m 2 ) 0 0 … 0 − ϕ ⁡ ( m 1 ) ​ ϕ ​ ( m 2 ) 2 ​ ϕ 2 ​ ( m 2 ) − ϕ ⁡ ( m 2 ) ​ ϕ ​ ( m 3 ) 0 … 0 0 − ϕ ⁡ ( m 2 ) ​ ϕ ​ ( m 3 ) 2 ​ ϕ 2 ​ ( m 3 ) − ϕ ⁡ ( m 3 ) ​ ϕ ​ ( m 4 ) … 0 ⋮ 0 0 0 0 … 2 ​ ϕ 2 ​ ( m n ) ) \begin{array}[]{l}\mathbf{V}^{-1}=(n+1)(n+2)\\ \times\left(\begin{array}[]{cccccc}2\phi^{2}\left(m_{1}\right)&-\phi\left(m_{1}\right)\phi\left(m_{2}\right)&0&0&\ldots&0\\ -\phi\left(m_{1}\right)\phi\left(m_{2}\right)&2\phi^{2}\left(m_{2}\right)&-\phi\left(m_{2}\right)\phi\left(m_{3}\right)&0&\ldots&0\\ 0&-\phi\left(m_{2}\right)\phi\left(m_{3}\right)&2\phi^{2}\left(m_{3}\right)&-\phi\left(m_{3}\right)\phi\left(m_{4}\right)&\ldots&0\\ \vdots&&&&&\\ 0&0&0&0&\ldots&2\phi^{2}\left(m_{n}\right)\end{array}\right)\end{array} (58)

## Appendix F Multivariate Statistics