ios - LAPACK sgelss 用法

Question

我想使用 LAPACK 框架中的 sgelss 函数来求解超定线性方程组。

minimize || A*x-b||

为此我选择使用 LAPACK 的 sgelss 函数，头文件可以在这里找到:SGELSS

我的问题是我似乎用错了，因为结果不合适。

一个例子，我使用这些矩阵:

double a[] =
    {
       25,5,1,
       16,4,1,
        9,3,1,
        4,2,1,
        1,1,1,
    };
double b[] =
     {
       31,21,13,7,3
     }

然后我会调用:

sgelss_(&(5), &(3), &(1), a, (5), b, &(5), 3, -1, &(3), malloc(3000), &(3000), &output);

[[ 因为 sgelss 期望大多数输入变量不是普通变量而是指针，所以我在这里使用 &(...) 来跳过实例化变量，然后为了这里的可读性而引用它们的内存位置! ]]

我希望之后 b 会是 [1,1,1,XXXXX] 因为我的输入 a 和 b 是这样设置的。不幸的是，情况并非如此。

我也尝试过旋转 a(切换行和列)，但没有成功。

Answer 1

您发布的代码中存在很多问题；我假设它是从你真正做的事情中抽象出来的，因为它显然永远无法按原样编译。有两个主要问题，我怀疑其中一个或两个是导致您的实际代码出现问题的原因:

首先，您尝试将 sgelss_ 与 double 数据一起使用。 sgelss_ 对float(单精度)进行操作。如果您有 double 据，则需要改用 dgelss_。

其次，您的参数描述了一个 5x3 矩阵，但请记住 LAPACK 对矩阵元素使用列优先排序。这意味着您的代码所描述的矩阵是:

25  1  2 
 5  9  1
 1  3  1
16  1  1
 4  4  1

不知怎的，我怀疑那不是你真正想要的矩阵。您更有可能需要矩阵:

25  5  1
16  4  1
 9  3  1
 4  2  1
 1  1  1

我继续整理了一个工作版本，假设这是您实际尝试使用的矩阵:

#include <Accelerate/Accelerate.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[]) {
    double a[] = /* column major storage order! */
    {
        25, 16,  9,  4,  1,
         5,  4,  3,  2,  1,
         1,  1,  1,  1,  1,
    };
    double b[] =
    {
        31,21,13,7,3
    };

    // Setup parameters
    __CLPK_integer m = 5;
    __CLPK_integer n = 3;
    __CLPK_integer nrhs = 1;
    __CLPK_integer lda = 5;
    __CLPK_integer ldb = 5;
    double *s = malloc(3 * sizeof*s);
    double rcond = -1.0f; // use machine precision
    __CLPK_integer rank;
    __CLPK_integer info;
    // Query correct worksize
    double worksize;
    __CLPK_integer lwork = -1;
    dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, &worksize, &lwork, &info);
    // Allocate workspace
    lwork = worksize;
    double *work = malloc(lwork * sizeof *work);
    // Do computation
    dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, work, &lwork, &info);
    // Free workspace
    free(work);
    // Print result vector
    for (int i=0; i<3; ++i)
        printf("%g\t", b[i]);
    printf("\n");
    return 0;
}