# MATLAB code for regularized Gaussian density-ratio experiments

This directory contains MATLAB scripts associated with **Regularized Variational and Spectral Log-Density-Ratio Estimation in the Gaussian Location Model**. The experiments use the whitened Gaussian location model

```text
q = N(0, I),        p = N(Delta, I),        s = ||Delta||^2,
alpha_p = m / n_p,  alpha_q = m / n_q.
```

The plotted quantities are excess risks under the paper's population criteria:

- `J`: log-normalized KL variational criterion.
- `K`: Fenchel/exponential KL criterion, which also checks the learned normalizer.
- `L`: two-potential spectral lower-bound criterion, used for the spectral estimator.

Several scripts separate numerical computation from plotting. Files named `*_simu.m` compute data and save `.mat` files. The corresponding plotting scripts load those files and write `.eps` figures.

## Archive status

This archive is **not a complete, directly runnable reproduction package**.

- Code is present for Figures 1, 2, 4, and 5.
- Simulation code is present for Figures 6, 7, and 8, but the plotting helper functions used by those figures are absent.
- No script is present for Figure 3, Figure 9, or Figure 10.
- The paper states that MATLAB code is available for all figures, so the archive is incomplete relative to that statement.
- The checks below are based on static inspection of the MATLAB files and visual/textual comparison with the supplied PDF. The scripts were not executed in this review.

## Figure-to-script map

| Paper figure | Script(s) | Output/data | Status and description |
|---|---|---|---|
| Figure 1 | `figure1_figure2.m` | `weaksignals.eps` | **Code present.** Fixed-signal first-order expansion of the unregularized risks. The three panels show the population term, the coefficient of `alpha_p`, and the coefficient of `alpha_q` for `J`, `K`, and `L`. |
| Figure 2 | `figure1_figure2.m` | `weaksignals_borders.eps` | **Code present, but an optional-looking template is required by an unguarded call.** Plots the zero-level sets of `Q_unreg`, `Q_J,opt`, and `Q_K,opt` in the weak-signal proportional regime. |
| Figure 3 | none | none | **Missing from archive.** Feasibility phase diagram for unregularized variational and spectral estimation at `s = 1`. |
| Figure 4, left | `figure4_simu.m`, then `figure4.m` | `fixed_regularization.mat`, `fixed_regularization.eps` | **Code present.** Near-unregularized comparison at `s = 1`, `tau = zeta = 1e-8`. Plots `log R_J^spec - log R_J^var` over approximately `[0,0.25]^2`, with the zero contour. |
| Figure 4, right | `figure4_simu_smalls.m`, then `figure4_smalls.m` | `fixed_regularization_smalls.mat`, `fixed_regularization_smalls.eps` | **Code present.** Same comparison at `s = 0.1`, `tau = zeta = 1e-8`, over approximately `[0,0.025]^2`. |
| Figure 5, top row | `figure5_simu.m`, then `figure5.m` | `best_regularization.mat`, `best_regularizationJ.eps`, `best_regularizationK.eps` | **Code present.** Optimized-regularization comparison at `s = 1` under `J` and `K`. |
| Figure 5, bottom row | `figure5_simu_smalls.m`, then `figure5_smalls.m` | `best_regularization_smalls.mat`, `best_regularizationJ_smalls.eps`, `best_regularizationK_smalls.eps` | **Code present.** Optimized-regularization comparison at `s = 0.1` under `J` and `K`. |
| Figure 6 | `figure6_simu.m`, then `figure6.m` | `sweep_m.mat`, intended `sweep_m.eps` | **Simulation code present; plotting dependency missing.** Dimension sweep at `s = 1`, `alpha_p = alpha_q = 0.125`, `tau = 1e-3`, `zeta = 1e-4`, with 32 replications. The intended plot compares empirical results with deterministic equivalents. |
| Figure 7 | `figure7_figure8_simu.m`, then `figure7_figure8.m` | `alpha_sweeps.mat`, intended `alphap_sweeps.eps` | **Simulation code present; plotting dependency missing.** Sweeps `alpha_p` over `[0,1]` with `alpha_q = 0.1`, using `s = 1`, `mVar = 1000`, `mSpec = 100`, 16 replications, `tau = 1e-2`, and `zeta = 1e-4`. |
| Figure 8 | `figure7_figure8_simu.m`, then `figure7_figure8.m` | `alpha_sweeps.mat`, intended `alphaq_sweeps.eps` | **Simulation code present; plotting dependency missing.** Sweeps `alpha_q` over `[0,1]` with `alpha_p = 0.05`, under the same settings as Figure 7. |
| Figure 9 | none | none | **Missing from archive.** Nuclear-penalty spectral estimator. The figure itself shows normalized criterion-`J` excess risk `2 R_J / s`: empirical performance with a shaded band versus the deterministic equivalent, as `alpha_p` is swept over `[0,1]` with `alpha_q = 0.1`. Appendix H.4 states `m = 100`. The prose below the figure incorrectly calls this an `alpha_q` sweep; the axis and caption show an `alpha_p` sweep. |
| Figure 10 | none | none | **Missing from archive.** Deterministic comparison of optimized nuclear-penalty and ridge spectral estimators at `s = 1`. The displayed quantity is `log R_nuclear - log R_ridge`; negative values favor the nuclear penalty, positive values favor ridge, and the black curve appears to be the zero contour. The paper does not identify which risk criterion is used. |

## Recommended run order

Run MATLAB from the directory containing these `.m` files.

```matlab
% Figures 1 and 2
figure1_figure2

% Figure 4
figure4_simu
figure4
figure4_simu_smalls
figure4_smalls

% Figure 5
figure5_simu
figure5
figure5_simu_smalls
figure5_smalls

% Figure 6: simulation works independently; plotting requires a missing helper
figure6_simu
figure6

% Figures 7 and 8: simulation works independently; plotting requires a missing helper
figure7_figure8_simu
figure7_figure8
```

Figures 3, 9, and 10 cannot be regenerated from this archive because their scripts are absent.

## Script details

### `figure1_figure2.m`

For Figure 1, the script evaluates fixed-signal first-order terms as functions of `s`. The variational expansion uses the affine closed form. The spectral expansion uses population deterministic formulas with 222-point Gauss-Legendre quadrature over `rho in [0,1]` and finite-difference derivatives with respect to `alpha_p` and `alpha_q`.

For Figure 2, the script plots the analytical zero contours of

```text
Q_unreg(beta_p,beta_q),
Q_J,opt(beta_p,beta_q),
Q_K,opt(beta_p,beta_q).
```

The call `openfig('singleplot.fig')` before Figure 2 is not protected by `try/catch`, while `singleplot.fig` is absent from the archive. Remove that call, guard it, or supply the template before running the complete script.

### `figure4_simu.m` and `figure4.m`

These produce the left panel of Figure 4. The computation uses

```text
s = 1,
tau = 1e-8,
zeta = 1e-8,
alpha_p in [0,0.25],
alpha_q in [0,0.25].
```

The code constructs the `alpha_q` grid by taking the union of a dense interval near zero and the `alpha_p` grid. Consequently, some comments and progress messages that state `alpha_q in [0,0.05]` are stale; the computed and plotted range extends to `0.25`.

### `figure4_simu_smalls.m` and `figure4_smalls.m`

These produce the right panel of Figure 4 with

```text
s = 0.1,
tau = 1e-8,
zeta = 1e-8,
alpha_p in [0,0.025],
alpha_q in [0,0.025].
```

As above, comments and progress messages that state `alpha_q in [0,0.005]` do not describe the final union grid.

### `figure5_simu.m` and `figure5.m`

These produce the top row of Figure 5 for `s = 1`. The computation covers `(alpha_p, alpha_q) in [0,0.5]^2` and searches

```text
regularization = 10^(-e),
e = -3, -2.875, ..., 6,
```

separately for each estimator and criterion. The spacing is `1/8` on the base-10 logarithmic exponent scale. The top-of-file comment saying `e = 1,...,6` is stale; the executable assignment is `-3:0.125:6`.

### `figure5_simu_smalls.m` and `figure5_smalls.m`

These repeat the Figure 5 calculation for `s = 0.1` and save `best_regularization_smalls.mat`.

### `figure6_simu.m` and `figure6.m`

The executable settings agree with the paper:

```text
s = 1,
alpha_p = alpha_q = 0.125,
tau = 1e-3,
zeta = 1e-4,
32 replications,
m = [50, 100:100:1000, 1200:200:2000, 2400, 2800].
```

The header comments in `figure6_simu.m` still state `alpha_p = alpha_q = 0.05` and `m = 100:100:1000`; those comments are stale. The simulation script saves the full replication arrays, means, and standard errors. The paper figure uses medians and interquartile bands, which are expected to be constructed by the missing plotting helper.

### `figure7_figure8_simu.m` and `figure7_figure8.m`

The executable settings agree with the paper:

```text
s = 1,
mVar = 1000,
mSpec = 100,
20 grid points,
16 replications,
tau = 1e-2,
zeta = 1e-4,
256 Gauss-Legendre nodes for the spectral deterministic equivalent.
```

The file headers incorrectly state `tau = 1e-4` and `zeta = 1e-5`, and the progress messages call the outputs “Figure 4” and “Figure 5.” The executable values and sweep definitions correspond to paper Figures 7 and 8.

## Supporting file

### `asymptotic_core_matlab.m`

Shared deterministic-equivalent dispatcher used by the Figure 4 and Figure 5 simulation scripts. It provides Gauss-Hermite and Gauss-Legendre quadrature, the variational entropy-dual solver, the spectral ridge-resolvent solver, metric conversion, and best-over-grid utilities.

The primary function is currently declared as

```matlab
function varargout = section7_asymptotic_core_matlab(action, varargin)
```

but the file is named `asymptotic_core_matlab.m` and the callers invoke `asymptotic_core_matlab(...)`. Align the primary function name with the filename, or rename the file and update every call. The current mismatch is not portable MATLAB code.

## Missing dependencies and files

The following referenced files are absent:

```text
plotMSweepQuantiles_outside.m    % required to plot Figure 6
plotSweepQuantiles_outside.m     % required to plot Figures 7 and 8
singleplot.fig                   % optional in most scripts, but unguarded for Figure 2
threeplots.fig                   % optional because its call is guarded
```

Several comments also refer to obsolete filenames such as `figure4_left_simu.m`, `figure4_right_simu.m`, and `best_regularization_contours_matlab.m`.

## Generated files

When all required dependencies are available, the included scripts are intended to generate

```text
weaksignals.eps
weaksignals_borders.eps
fixed_regularization.mat
fixed_regularization.eps
fixed_regularization_smalls.mat
fixed_regularization_smalls.eps
best_regularization.mat
best_regularizationJ.eps
best_regularizationK.eps
best_regularization_smalls.mat
best_regularizationJ_smalls.eps
best_regularizationK_smalls.eps
sweep_m.mat
sweep_m.eps
alpha_sweeps.mat
alphap_sweeps.eps
alphaq_sweeps.eps
```

No output filenames or generating scripts for Figures 3, 9, or 10 are present.

## Discrepancies found

1. **Figure 9 variable mismatch in the paper.** The figure axis and caption specify an `alpha_p` sweep with `alpha_q = 0.1`; Appendix H.4 prose says it is an `alpha_q` sweep. The figure and caption are mutually consistent, so `alpha_p` is the more likely intended variable. Confidence: approximately 99%.
2. **Figure 10 does not identify the criterion.** The plot and caption use a generic `R`, not `R_J`, `R_K`, or `R_L`. With no Figure 10 script, the criterion and optimization grid cannot be checked. Confidence that this is genuinely unspecified: approximately 98%.
3. **Figure 10's “wide range” claim is asymmetric.** The plotted `alpha_q` range is broad (`0` to `0.5`), but the `alpha_p` axis is only approximately `0` to `5e-4` because of the `x 10^-3` scale. Describing the plot as covering a wide range of both aspect ratios is misleading. Confidence: approximately 95%.
4. **Figure 4 body text is incomplete for the right panel.** The surrounding text describes `s = 1` and `[0,1/4]^2`, while the caption, plot, and code also include the `s = 0.1` panel on `[0,0.025]^2`. Confidence: approximately 99%.
5. **The archive is incomplete relative to the paper's reproducibility statement.** Scripts for Figures 3, 9, and 10 and plotting helpers for Figures 6-8 are absent. Confidence: 100% from the supplied archive inventory.
6. **Several code comments do not match executable parameters.** This affects the Figure 4 grid descriptions, Figure 5 regularization range, Figure 6 aspect ratios and dimension list, and Figure 7-8 ridge values and figure numbers. The executable values generally match the paper more closely than the comments. Confidence: approximately 99%.

## Notes

- Run scripts from the directory containing the `.m` files so MATLAB can find local files and saved `.mat` data.
- Plotting scripts overwrite their `.eps` outputs.
- Simulation scripts overwrite their `.mat` outputs.
- Figures 6-8 cannot be plotted without the missing helper functions.
- Figures 3, 9, and 10 cannot be reproduced from this archive as supplied.



This package is provided free for non-commercial use under the terms of the GNU General Public License. If you have any questions or comments regarding this package, or if you want to report any bugs, please send an e-mail to francis.bach@inria.fr. 


