c - 求解使 5x5 对称矩阵奇异所需的值

Question

我有一个看起来像这样的 5x5 矩阵(具体来说，它是嵌入到平面中的 4 点形状的边长 Cayley-Menger 矩阵):

如果所有其他值已知，我正在尝试求解使矩阵奇异的值(因此行列式等于零)，即:

a: 2.357086
b: 2.017018
c: 2.681643
d: <unknown>
e: 2.900140
f: 3.261168

在上面的例子中，我只是出于演示目的而省略了 d(不过，如果要求解矩阵，答案将是 1.597908)。可能会缺少 a、b、c、d、e 或 f，但只有一个，而且永远只有一个。

我如何以相当通用的方式在代码中解决这个问题？最终，我希望创建一个函数，在其中可以传入 5 个必需变量和一个索引参数，该索引参数确定要求解的未知数，然后返回所有可能的答案(可能是 0、1 或 2 个可能的答案) - 像这样的东西:

vector2 SolveForMissingSide(float sides[], int index)
{
    // Do something here, return 0, 1, or 2 answers
}

float sides[] = {2.357086, 2.017018, 2.681643, 0, 2.900140, 3.261168};
//                                             ^
//                        4th side set to zero here because it's unknown

// Call the solver to solve for the 4th unknown side, or "d" in this case

vector2 missingSides = SolveForMissingSides(sides, 4);

printf("The length of side d is either: %f or %f\n", missingSides[0], missingSides[1]);

不幸的是，我在这里无法访问 Eigen 或 numpy 之类的东西，因为我正在尝试为嵌入式设备编写一个函数，该函数对类似 C 的设施的访问相当有限(常见的数学函数，如 sqrt、sin/cos/tan 等都可用)。

Answer 1

以下答案一般来说不太相关， 但由于提出了一个具有挑战性的问题，所以就这样了。

我首先在 Mathematica 中计算行列式:

det = {{0, a^2, e^2, d^2, 1},{a^2, 0, b^2, f^2,1}, {e^2, b^2, 0, c^2, 1}, {d^2, f^2, c^2, 0, 1}, {1, 1, 1, 1, 0}} //Det

经过一些操作，C代码生成并使用 a2 = a^2, b2 = b^2, .... 行列式简化为:

det = -2*(a2*a2*c2 - a2*b2*c2 - a2*b2*d2 + a2*b2*e2  + a2*c2*c2 - a2*c2*d2 - a2*c2*e2 - a2*c2*f2 +
 a2*d2*f2 - a2*e2*f2 + b2*b2*d2 - b2*c2*d2 + b2*c2*f2 + b2*d2*d2 - b2*d2*e2 - b2*d2*f2 -
 b2*e2*f2 + c2*d2*e2 - c2*e2*f2 - d2*e2*f2 + e2*e2*f2 + e2*f2*f2)

由于我们需要求解方程 det = 0，因此我们忽略 -2 乘法因子，剩下的是 22 个项由 +1 或 -1 因子组成，并且 六个变量中的三个相乘 a2, b2, ...。

这些术语可能存在逻辑，因此可以以智能方式生成它们，例如与变量排列的签名相关的东西。不过，我并没有花费太多精力去寻找这种逻辑，而是只是从 Mathematica 中对它们进行了硬编码。

因此，这些项被实现为数组 int terms[22][4]，每个项对应的四个整数是:首先是 +1 或 -1系数，然后是变量的三个索引 按照a2->0、b2->1、...f2->5的顺序。

由于我们知道六个变量中的五个a2, b2, ..., f2， 很容易找到二阶的三个系数 det = 0 的等式。

这是实现 - 抱歉使用了 double 而不是 float 以及不同的函数签名(但是，有四种可能的解决方案，我不知道哪一种选择一项)， 但我确信OP可以适应它:

#include <stdio.h>
#include <math.h>

const int terms[22][4] = {
        {+1, 0, 0, 2}, //+a2*a2*c2
        {-1, 0, 1, 2}, //-a2*b2*c2
        {-1, 0, 1, 3}, //-a2*b2*d2
        {+1, 0, 1, 4}, //+a2*b2*e2
        {+1, 0, 2, 2}, //+a2*c2*c2
        {-1, 0, 2, 3}, //-a2*c2*d2
        {-1, 0, 2, 4}, //-a2*c2*e2
        {-1, 0, 2, 5}, //-a2*c2*f2
        {+1, 0, 3, 5}, //+a2*d2*f2
        {-1, 0, 4, 5}, //-a2*e2*f2

        {+1, 1, 1, 3}, //+b2*b2*d2
        {-1, 1, 2, 3}, //-b2*c2*d2
        {+1, 1, 2, 5}, //+b2*c2*f2 
        {+1, 1, 3, 3}, //+b2*d2*d2 
        {-1, 1, 3, 4}, //-b2*d2*e2 
        {-1, 1, 3, 5}, //-b2*d2*f2 
        {-1, 1, 4, 5}, //-b2*e2*f2 

        {+1, 2, 3, 4}, //+c2*d2*e2 
        {-1, 2, 4, 5}, //-c2*e2*f2 

        {-1, 3, 4, 5}, //-d2*e2*f2 

        {+1, 4, 4, 5}, //+e2*e2*f2
        {+1, 4, 5, 5}, //+e2*f2*f2
};

void solveDet(int unknownIdx, double knownValues[5], int *nSolutions, double solutions[4]){
    // all known vales squared, and 1 for the unknown variable
    double allSqrValues[6];
    for(int i = 0; i < 6; i++){
        allSqrValues[i] = i < unknownIdx ? knownValues[i] :
                i == unknownIdx ? 1 : knownValues[i-1];
        allSqrValues[i] *= allSqrValues[i];
    }
    double a = 0, b = 0, c = 0; // coeffs of second degree equation
    for(int iTerm = 0; iTerm < 22; iTerm++){
        const int* term = terms[iTerm];
        double term_value = term[0] * allSqrValues[term[1]] *
                allSqrValues[term[2]] * allSqrValues[term[3]];
        // the number of occurrences of the unknown variable in the term
        int nOccurrences = (term[1] == unknownIdx) +
                (term[2] == unknownIdx) + (term[3] == unknownIdx);
        if(nOccurrences == 2){ //x^2
            a += term_value;
        }
        else if(nOccurrences == 1){ //x^1
            b += term_value;
        }
        else{ //x^0 - no x
            c += term_value;
        }
    }
    // The solution of the 2nd degree equation
    double sqrtDelta = sqrt(b*b - 4*a*c),
        x1 = (-b - sqrtDelta)/2/a,
        x2 = (-b + sqrtDelta)/2/a;
    // return ± the square root of the positive solutions - 0, 2, or 4 solutions
    *nSolutions = 0;
    if(x1 >= 0){
        solutions[(*nSolutions)++] = -sqrt(x1);
        solutions[(*nSolutions)++] = sqrt(x1);
    }
    if(x2 >= 0){
        solutions[(*nSolutions)++] = -sqrt(x2);
        solutions[(*nSolutions)++] = sqrt(x2);
    }

    //  // verify det = 0
    //  allSqrValues[unknownIdx] = x2;
    //  double det = 0;
    //  for(int iTerm = 0; iTerm < 22; iTerm++){
    //      const int *term = terms[iTerm];
    //      //printf("%2d %d %d %d\n", term[0], term[1], term[2], term[3]);
    //      det += term[0] * allSqrValues[term[1]] *
    //             allSqrValues[term[2]] * allSqrValues[term[3]];
    //  }
    //  printf("det=%g\n", det);
};


int main(){
    //a: 2.357086
    //b: 2.017018
    //c: 2.681643
    //d: <unknown>
    //e: 2.900140
    //f: 3.261168
    double knownValues[5] = {2.357086, 2.017018,
                             2.681643, 2.900140, 3.261168};
    int unknownIndex = 3;
    int nSolutions;
    double solutions[4];
    solveDet(unknownIndex, knownValues, &nSolutions, solutions);
    if(nSolutions == 0){
        printf("No solutions\n");
    }
    else{
        printf("Solutions: ");
        for(int iSol = 0; iSol < nSolutions; iSol++){
            printf("%g ", solutions[iSol]);
        }
        printf("\n");
    }
    return 0;
}

Online GDB link