计算科学 - 当 Inverse 不存在但行列式提供非零结果时克服浮点问题 - 吾爱随笔录

我有一个表达式（假设矩阵 A 的行列式）以符号形式表示，用 2 个决策变量 x1、x2 和 2 个参数 q1 和 q2 表示。对于不同的起点和参数值组合，我使用 fmincon 将其最小化。

对于 x1、x2、q1、q2 的某些组合，矩阵 A 是不可逆的，所以我希望行列式评估为零，但 Matlab 只给出一个非常小的数字，如 9.7166e-91。

如果我用不存在逆的数值创建了矩阵 A，然后取逆，它确实会给我通常的警告“警告：矩阵接近奇异或缩放错误。结果可能不准确。” 这个警告非常好，因为我知道可以丢弃输出。

但我需要使用符号表达式，因为这是需要提供给 fmincon 的。

那么我应该使用什么方法来知道应该丢弃符号表达式的数值评估的某些输出？就像很高兴在不存在逆时拥有“NaN”一样。希望有一个不会减慢我的程序的技巧

请注意，数字设置为 32，所以我应该使用它并推断如果行列式结果介于 -1E-32 和 1E-32 之间，那么结果应该被丢弃吗？（这种方法似乎不正确）

欢迎任何建议。

%Providing code below in order to replicate easily;

syms x1 x2 q1 q2;% x1 and x2 are the decision variables while q1 and q2 are the parameters. For different values of parameters I'm trying to minimize a function using fmincon


%matrix A is DesignMatrixMult. Example provided below for illustrative purposes but its form can change depending on user input;
 DesignMatrixMult = [                    ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 - (exp(-q1*x1) - exp(-q2*x1))/(q1 - q2) + (q1*x1*exp(-q1*x1))/(q1 - q2))^2 + ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 - (exp(-q1*x2) - exp(-q2*x2))/(q1 - q2) + (q1*x2*exp(-q1*x2))/(q1 - q2))^2, - ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 + (q1*x1*exp(-q2*x1))/(q1 - q2))*((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 - (exp(-q1*x1) - exp(-q2*x1))/(q1 - q2) + (q1*x1*exp(-q1*x1))/(q1 - q2)) - ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1
    - q2)^2 + (q1*x2*exp(-q2*x2))/(q1 - q2))*((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 - (exp(-q1*x2) - exp(-q2*x2))/(q1 - q2) + (q1*x2*exp(-q1*x2))/(q1 - q2));
     - ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 + (q1*x1*exp(-q2*x1))/(q1 - q2))*((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 - (exp(-q1*x1) - exp(-q2*x1))/(q1 - q2) + (q1*x1*exp(-q1*x1))/(q1 - q2)) - ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1
    - q2)^2 + (q1*x2*exp(-q2*x2))/(q1 - q2))*((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 - (exp(-q1*x2) - exp(-q2*x2))/(q1 - q2) + (q1*x2*exp(-q1*x2))/(q1 - q2)),                                        ((q1*(exp(-q1*x1) - exp(-q2*x1)))/(q1 - q2)^2 + (q1*x1*exp(-q2*x1))/(q1 - q2))^2 + ((q1*(exp(-q1*x2) - exp(-q2*x2)))/(q1 - q2)^2 + (q1*x2*exp(-q2*x2))/(q1 - q2))^2];


% determinant is calculated below;  
DCritn = simplify(simplify(det(DesignMatrixMult),'IgnoreAnalyticConstraints', true,'steps',500),'full');
%creating an anonymous function handle to DCritn;
DCritnF = matlabFunction(DCritn,'vars',{x1,x2,q1,q2},'file','');



% few lines of code for generating different possible values of the parameter q1 and q2;
% As an example, one set of values is provided below;
q1=0.9287; q2 = 0.83;

对于上述参数值，函数 DCritnF 使用 fmincon 最小化；

我面临的挑战是，对于 x1 和 x2 的某些组合，即使我知道 DesignMatrixMult（矩阵 A）的逆不存在，DCritnF 也会给出非零结果。例如，使用 x1 = 42.5 和 x2 = 95 的值调用 DCritnF，它返回 9.7166e-91

另一方面，如果我执行以下操作，

TrialCheck = double(subs(DesignMatrixMult));
TrialCheckInverse = inv(TrialCheck);

警告：矩阵接近奇异或严重缩放。结果可能不准确。RCOND = 5.404300e-18。答案=

1.0e+42 *

5.4788 -1.3938

-1.3938 0.3546

还有一个小问题；如果我做了

TrialCheck_Alternate = double(subs(DesignMatrixMult))
TrialCheckDeterm_Alternate = det(TrialCheck_Alternate)

然后我没有收到警告？如果对于矩阵的逆它会发出警告，那么它不应该在找到该矩阵的行列式时也发出警告吗？

TrialCheck_Alternate =

1.0e-25 *

0.0139 0.0546 0.0546 0.2146

TrialCheckDeterm_Alternate =

3.9175e-69