Reproduce Auto-modeling on Many-Normal-Means
\[\newcommand\bbP}\notag\]Lijun Wang (2022.07.31)
Here is my attempt on reproducing the results of auto-modeling on the many-Normal-means example.
The main algorithm is
- Take another perspective of (3.3), it is just to minimize $G_{\hat\bbP}(\theta, \lambda)$. Thus, we can directly take advantage with existing nonlinear programming solvers, such as Ipopt (Interior Point Optimizer).
- For (3.4) and (3.5), since $\theta_2,\cdots,\theta_m\ge 0$, then
and hence
\[\frac{\partial \pi^{(n)}(\theta,\lambda)}{\partial \theta} = % \begin{bmatrix} % 0 \\ % \lambda_2\\ % \vdots\\ % \lambda_m\\ % 0\\ % \vdots\\ % 0 % \end{bmatrix} \begin{bmatrix} 0 & \lambda_2 & \lambda_3 & \cdots & \lambda_m & 0 & \cdots & 0 \end{bmatrix}^T\]And note that $\lambda_i\ge 0$, it follows that \(\lambda_i = \max(0, d_i), i=2,\ldots,m\) where $d_i$ is the $i$-th element of \(E_{(X, Y)\sim \bbP}[\frac{\partial L(\theta\mid X, Y)}{\partial \theta}] - \frac 1n\sum_{i=1}^n\frac{\partial L(\theta\mid x_i, y_i)}{\partial \theta}\) For the many-Normal-means example,
\[\begin{align*} L_i \triangleq L(\theta\mid y_i) &= -\log \sum_{k=1}^m \alpha_k \exp\left(-\frac{(y_i-\eta_k)^2}{2}\right)\\ &= -\log \sum_{k=1}^m \theta_{m+k} \exp\left(-\frac{(y_i-\sum_{j=1}^k\theta_j)^2}{2}\right)\,, \end{align*}\]then
\(\frac{\partial L_i}{\partial \theta_\ell} = \begin{cases} -\frac{\sum_{k=\ell}^m\theta_{m+k}\exp\left(-\frac{(y_i-\sum_1^k\theta_j)^2}{2}\right)(y_i-\sum_1^k\theta_j)}{\sum_{k=1}^m\theta_{m+k}\exp\left(-\frac{(y_i-\sum_1^k\theta_j)^2}{2}\right)} & \ell\le m\\ -\frac{\exp\left(-\frac{(y_i-\sum_1^{\ell-m}\theta_j)^2}{2}\right)}{\sum_{k=1}^{m}\theta_{m+k}\exp\left(-\frac{(y_i-\sum_1^k\theta_j)^2}{2}\right)}& \ell\ge m+1 \end{cases}\) For simplicity, I skip the step of approximation of $\bbP$ by bootstrap samples, and instead directly use the true $\bbP$.
My pseudo Julia code is as follows,
function auto_modeling()
for i = 1:N
θold = θ
λold = λ
λ = sol_λ_given_θ(y, θ, ...)
θ = sol_θ_given_λ(y, λ, ...)
if (‖ θold - θ ‖ < tol) and (‖ λold - λ ‖ < tol)
break
end
end
end
function sol_λ_given_θ(y, θ, ...)
# Equations (3) (4) (5)
end
function sol_θ_given_λ(y, λ, ...)
model = Model(Ipopt.Optimizer)
# express the optimization problem `min G` in language of JuMP and Ipopt
optimize!(model)
end
As for g-modeling, call deconv function from the R package deconvolveR.
Like in the paper, I repeat 200 times, then report the average mean square error. The results are as follows:
Compared to the results in the paper,
- Except for the results of g-modeling, others are close to the reported results in the paper.
- In my experiments, g-modeling can outperform auto-modeling, but I did not deliberately select its parameters. The paper also did not discuss how they chose the parameters for this method. So g-modeling might be better than the reported performance in the paper.
- Currently, the program is relatively slow, and the computational burden is mainly Step 3 in the coordinate descent algorithm. Are there any speed-up strategies?