Numerical error with `MPR`, `t = 0.01` and `m = 20`

[AlexisDerumigny]

With the attached dataset dataset_numerical_error.csv, we obtain an error since the quantity h that is computed here

UniversalShrink/R/ridge_higher_order_shrinkage.R

Lines 450 to 453 in 7f23ab0

	h[j + 1] = h_hat_jp1_t(h_hat_until_j = h[1:j],
	v_hat_until_j = v[1:j],
	j = j,
	Bell_polynomials = Bell_polynomials)

is not finite.

This happens when running the following code:

X = read.csv("dataset_numerical_error.csv", sep = ";", dec = ".", header  = FALSE) |> 
  as.matrix()

estimator = UniversalShrink::MPR_higher_order_shrinkage(
    X = X, t = 0.01, m = 20, verbose = 3)

Here are some more detailed debugging information.

Starting `MPR_higher_order_shrinkage`...
*  n =  100 
*  p =  200 
*  centered case
*  c_n =  2.020202 

*  t =  0.01 
*  m =  20 
Starting values: 
*  q1 =  5.434845 
*  q2 =  45.41597 

Estimation of v(t0) =  0.2477956 

Estimation of the derivatives of v:
 [1]  -1.448527e-01   2.557478e-01  -8.150259e-01   3.812940e+00  -2.827592e+01
 [6]   3.579545e+03  -2.718545e+06   2.444940e+09  -2.444231e+12   2.690455e+15
[11]  -3.241584e+18   4.206028e+21  -5.892296e+24   8.823948e+27  -1.413418e+31
[16]   2.402487e+34  -4.326465e+37   8.209085e+40  -1.642912e+44   3.452907e+47
[21]  -7.594814e+50   1.745789e+54  -4.186255e+57   1.047818e+61  -2.724067e+64
[26]   7.353747e+67  -2.057596e+71   5.966502e+74  -1.789339e+78   5.549470e+81
[31]  -1.777162e+85   5.862187e+88  -1.994182e+92   6.974020e+95  -2.512183e+99
[36]  9.288270e+102 -3.532097e+106  1.376800e+110 -5.505985e+113  2.258029e+117
[41] -9.481565e+120  4.080574e+124 -1.795249e+128  8.076528e+131 -3.714669e+135
[46]  1.746462e+139 -8.381833e+142  4.106897e+146 -2.053038e+150  1.046954e+154
[51] -5.445669e+157  2.887760e+161 -1.558238e+165  8.572075e+168 -4.799304e+172
[56]  2.736208e+176 -1.586612e+180  9.364709e+183 -5.616417e+187  3.428211e+191
[61] -2.124661e+195  1.338775e+199 -8.565044e+202  5.566158e+206 -3.675160e+210
[66]  2.462523e+214 -1.674267e+218  1.155062e+222 -8.088074e+225  5.739871e+229
[71] -4.134112e+233  3.018363e+237 -2.232570e+241  1.674509e+245 -1.272967e+249
[76]  9.803544e+252 -7.644955e+256  6.039782e+260 -4.832344e+264

Estimation of h:
 [1]             NA   6.903563e+00   4.207284e+01   2.042733e+02   9.026184e+02
 [6]  -3.698127e+02   3.362859e+06  -2.613274e+09   2.011618e+12  -1.528931e+15
[11]   1.151661e+18  -8.634341e+20   6.384071e+23  -4.708963e+26   3.445512e+29
[16]  -2.518594e+32   1.830434e+35  -1.326906e+38   9.567664e+40  -6.898498e+43
[21]   4.960342e+46  -3.554285e+49   2.540925e+52  -1.812532e+55   1.294084e+58
[26]  -9.204755e+60   6.539636e+63  -4.636522e+66   3.286169e+69  -2.325078e+72
[31]   1.645127e+75  -1.162938e+78   8.199487e+80  -5.785682e+83   4.069917e+86
[36]  -2.867843e+89   2.014167e+92  -1.417392e+95   9.942607e+97 -6.977654e+100
[41]  4.895895e+103 -3.430115e+106  2.406286e+109 -1.683687e+112  1.177979e+115
[46] -8.238719e+117  5.763060e+120 -4.025552e+123  2.812037e+126 -1.962813e+129
[51]  1.369903e+132 -9.561859e+134  6.672848e+137 -4.644924e+140  3.240315e+143
[56] -2.255996e+146  1.572178e+149 -1.093936e+152  7.621341e+154 -5.297320e+157
[61]  3.690904e+160 -2.563382e+163  1.783807e+166 -1.238780e+169  8.609430e+171
[66] -5.986512e+174  4.156941e+177 -2.886239e+180  2.003693e+183 -1.391914e+186
[71]  9.647914e+188 -6.703295e+191  4.649010e+194 -3.221795e+197  2.235846e+200
[76] -1.550687e+203           -Inf            NaN            NaN            NaN

The problem is simply that the numbers become too big to represent.
One should probably switch to Rmpfr in those cases...

Similar errors can also be reproduced with smaller t and m, for example

estimator = UniversalShrink::MPR_higher_order_shrinkage(
    data_strange, t = 0.0001, m = 14, verbose = 3)