Utilities

hermite(x::Real, n::Int, probabilists::Bool=true)

Compute the Hermite polynomials of degree n. Compute the Probabilists' version by default.

The two definitions of the Hermite polynomials are each a rescaling of the other. Let $Heₙ(x)$ denote the Probabilists' version, and $Hₙ(x)$ the Physicists'. Then

\[H_{n}(x) = 2^{\frac{n}{2}} He_{n}\left(\sqrt{2} x\right)\]

\[He_{n}(x) = 2^{-\frac{n}{2}} H_{n}\left(\frac{x}{\sqrt{2}}\right)\]

normal_to_margin(d::UnivariateDistribution, x::Float64)

Convert samples from a standard normal distribution to a given marginal distribution.

_randn(n::Int, d::Int)

Fast parallel generation of standard normal samples.

_rmvn(n::Int, ρ::Matrix{Float64})

Fast parallel generation of multivariate standard normal samples.

get_coefs(margin::UnivariateDistribution, n::Int)

Get the $n^{th}$ degree Hermite Polynomial expansion coefficients for $F^{-1}[Φ(⋅)]$ where $F^{-1}$ is the inverse CDF of a probability distribution and Φ(⋅) is the CDF of a standard normal distribution.

Notes

The paper describes using Guass-Hermite quadrature using the Probabilists' version of the Hermite polynomials, while the package FastGaussQuadrature.jl uses the Physicists' version. Because of this, we need to do a rescaling of the input and the output:

\[\frac{1}{k!}\sum_{s=1}^{m}w_s H_k (t_s) F_{i}^{-1}\left[\Phi(t_s)\right] ⟹ \frac{1}{\sqrt{\pi} \cdot k!}\sum_{s=1}^{m}w_s H_k (t_s\sqrt{2}) F_{i}^{-1}\left[\Phi(t_s)\right]\]

Hϕ(x::Real, n::Int)

We need to account for when x is ±∞ otherwise Julia will return NaN for 0×∞

Gn0d(n::Int, A, B, α, β, σAσB_inv)

Calculate the $n^{th}$ derivative of G at 0 where $ρ_x = G(ρ_z)$

We are essentially calculating a double integral over a rectangular region

\[\int_{α_{r-1}}^{α_r} \int_{β_{s-1}}^{β_s} Φ(z_i, z_j, ρ_z) dz_i dz_j\]

(α[r], β[s+1]) +----------+ (α[r+1], β[s+1])
               |          |
               |          |
               |          |
  (α[r], β[s]) +----------+ (α[r+1], β[s])

Gn0m(n::Int, A, α, dB, σAσB_inv)

Calculate the $n^{th}$ derivative of G at 0 where $ρ_x = G(ρ_z)$

solve_poly_pm_one(coef)

Solve a polynomial equation on the interval [-1, 1].

npd_gradient(y::Vector{Float64}, λ₀::Vector{Float64}, P::Matrix{Float64}, b₀::Vector{Float64}, n::Int)

Return f(yₖ) and ∇f(yₖ) where

\[f(y) = \frac{1}{2} \Vert (A + diag(y))_+ \Vert_{F}^{2} - e^{T}y\]

and

\[\nabla f(y) = Diag((A + diag(y))_+) - e\]

npd_pca(b::Vector{Float64}, X::Matrix{Float64}, λ::Vector{Float64}, P::Matrix{Float64}, n::Int)

npd_pre_cg(b::Vector{Float64}, c::Vector{Float64}, Ω₀::Matrix{Float64}, P::Matrix{Float64}, ϵ::Float64, N::Int, n::Int)

Preconditioned conjugate gradient method to solve Vₖdₖ = -∇f(yₖ)

npd_precond_matrix(Ω₀::Matrix{Float64}, P::Matrix{Float64}, n::Int)

Create the precondition matrix used in solving the linear system Vₖdₖ = -∇f(yₖ) in the conjugate gradient method.

npd_set_omega(λ::Vector{Float64}, n::Int)

Used in creating the precondition matrix.

npd_jacobian(x::Vector{Float64}, Ω₀::Matrix{Float64}, P::Matrix{Float64}, n::Int)

Create the Generalized Jacobian matrix for the Newton direction step.

fast_pca!(X::Matrix{T}, λ::Vector{T}, P::Matrix{T}, n::Int) where T<:AbstractFloat