在阅读其他答案并进行调查后，我明白我的错误是浮点错误。对于我的用例，我还没有看到任何提供良好解决方案的答案。
上下文：在 r 中使用sSDR包中的matpower函数来帮助我标准化预测因子X数据。
我的数据X形状如下(实际数据集有更多的行和列)：
Press_mm_hg RH_out Windspeed Visibility Tdewpoint rv1 rv2
1 733.5 92 7.000000 63.00000 5.3 13.27543 13.27543
2 733.6 92 6.666667 59.16667 5.2 18.60619 18.60619
3 733.7 92 6.333333 55.33333 5.1 28.64267 28.64267
4 733.8 92 6.000000 51.50000 5.0 45.41039 45.41039
5 733.9 92 5.666667 47.66667 4.9 10.08410 10.08410
6 734.0 92 5.333333 43.83333 4.8 44.91948 44.91948
我想把var(X)提升到-0.5的力量.我这样做如下(从sSDR包matpower函数中提取代码)：
# options(scipen=999) helps you visualise the problem better.
X_var = var(X) # compute var of X.
X_var_sum = (var_X + t(var_X))/2 # Sum the result
tmp = eigen(X_var_sum) # Eigenvalue decomposition
e_values <- tmp$values # Get the values
e_vectors <- tmp$vectors # Get the vectors
e_vectors %*% diag(e_values^alpha) %*% t(e_vectors) # Raise the matrix to the power.
那么错误是什么？
错误发生在特征值分解期间。我的数据结构的某些内容导致特征值为负数。例如
> e_values
[1] 497.855846350015326606808230280876159667969
[2] 37.044927498529837350815796526148915290833
[3] 0.000000000000034235240904804834595748182
[4] 0.000000000000000006398983555805967132521
[5] 0.000000000000000000000000000000000385186
[6] -0.000000000000000003469446977025200723022
[7] -0.000000000000012918418921653726177566030
这对将特征值提高到 -0.5 的幂有进一步的影响。由于负数e_values^alpha具有NaN值。
可能的解决方案
将负特征值设置为 0。
0^(-0.5)是无限的。- 所以不可行，因为未来的计算依赖于此。
设置使值为正。(我当前的解决方案)
- 由于数字太小，由于错误，
abs(e_values)返回巨大的数字。
- 由于数字太小，由于错误，
将值设置为几乎 0，但刚好高于 0。
- 我认为这将完全扭曲数据
问题
有没有更好的方法来处理这个问题，或者可能有不同的特征值分解算法，它不太容易出现这些错误？
dput()X数据
`structure(list(Press_mm_hg = c(733.5, 733.6, 733.7, 733.8, 733.9,
734), RH_out = c(92, 92, 92, 92, 92, 92), Windspeed = c(7, 6.66666666666667,
6.33333333333333, 6, 5.66666666666667, 5.33333333333333), Visibility = c(63,
59.1666666666667, 55.3333333333333, 51.5, 47.6666666666667, 43.8333333333333
), Tdewpoint = c(5.3, 5.2, 5.1, 5, 4.9, 4.8), rv1 = c(13.275433157105,
18.6061949818395, 28.6426681675948, 45.4103894997388, 10.0840965518728,
44.9194842483848), rv2 = c(13.275433157105, 18.6061949818395,
28.6426681675948, 45.4103894997388, 10.0840965518728, 44.9194842483848
)), row.names = c(NA, 6L), class = "data.frame")
您可以通过结合使用zapsmall和扰动来获得输出.Machine$double.eps，但我会对输出持谨慎态度。
e_vectors %*% diag((zapsmall(e_values)+.Machine$double.eps)^alpha) %*% t(e_vectors)
[,1] [,2] [,3] [,4] [,5]
[1,] 6.706360e+07 -1.280137e-09 1.508856e+05 1.735184e+06 4.526567e+04
[2,] -1.280137e-09 6.710886e+07 2.793968e-09 5.913898e-08 -2.328306e-10
[3,] 1.508856e+05 2.793968e-09 6.660591e+07 -5.783946e+06 -1.508856e+05
[4,] 1.735184e+06 5.913898e-08 -5.783946e+06 5.934833e+05 -1.735184e+06
[5,] 4.526567e+04 -2.328306e-10 -1.508856e+05 -1.735184e+06 6.706360e+07
[6,] -3.874249e-04 5.444576e-09 1.291434e-03 1.485146e-02 3.874244e-04
[7,] -3.874249e-04 -5.029556e-09 1.291425e-03 1.485146e-02 3.874330e-04
[,6] [,7]
[1,] -3.874249e-04 -3.874249e-04
[2,] 5.444576e-09 -5.029556e-09
[3,] 1.291434e-03 1.291425e-03
[4,] 1.485146e-02 1.485146e-02
[5,] 3.874244e-04 3.874330e-04
[6,] 3.355443e+07 -3.355443e+07
[7,] -3.355443e+07 3.355443e+07