Amoiensis / Matrix_hub

A lib of Matrix operation for C language. (矩阵运算库--C语言)
Apache License 2.0
234 stars 53 forks source link

householder函数向量模为0情况下计算错误以及是否需要添加0 == _x->data[0]的情况 #11

Open 645770225 opened 1 year ago

645770225 commented 1 year ago

作者你好,感谢开源 householder函数使用过程中遇到一些问题,当待变换的向量为0向量时,会计算错误,我将householder函数前4行修改为:

    Matrix *H = NULL;
    Matrix *y = M_Zeros(_x->row, _x->column);
    MATRIX_TYPE x_norm = M_norm(_x, 2);
    if (x_norm < 1e-5f)
    {
        M_free(y);
        H = M_I(_x->row);
        return H;
    }
    y->data[0] = x_norm;
    Matrix *w = NULL;

当向量为0向量时,将H阵设为单位阵,计算就不会出现错误了, 此外我看参考文档https://max.book118.com/html/2017/0527/109650252.shtm 中提到sgn(a)在a=0的时候为0,所以将if分支修改为:

  if (_x->data[0] > 1e-5f)
    {
        w = M_add_sub(1, _x, -1, y);
        Matrix *temp_w = w;
        w = M_numul(w, 1 / M_norm(w, 2));
        M_free(temp_w);
    }
    else if (fabs(_x->data[0]) <= 1e-5f)
    {
        w = _x;
        w = M_numul(w, 1 / M_norm(w, 2));
    }
    else
    {
        w = M_add_sub(1, _x, 1, y);
        Matrix *temp_w = w;
        w = M_numul(w, 1 / M_norm(w, 2));
        M_free(temp_w);
    }

从而和文档上对应。

yukun1107 commented 7 months ago

作者你好,这部分修改代码有一个BUG需要修正一下。 else if (fabs(_x->data[0]) <= 1e-5f) { w = _x; w = M_numul(w, 1 / M_norm(w, 2)); } 中的w=_x应该改为copy,不然_x会在这个函数中提前被free。 应该改成: else if (fabs(_x->data[0]) <= 1e-5f) { w = Matrix_copy(_x); w = M_numul(w, 1 / M_norm(w, 2)); }

645770225 commented 7 months ago

作者你好,这部分修改代码有一个BUG需要修正一下。 else if (fabs(_x->data[0]) <= 1e-5f) { w = _x; w = M_numul(w, 1 / M_norm(w, 2)); } 中的w=_x应该改为copy,不然_x会在这个函数中提前被free。 应该改成: else if (fabs(_x->data[0]) <= 1e-5f) { w = Matrix_copy(_x); w = M_numul(w, 1 / M_norm(w, 2)); }

这里不会提前释放的,w在w = M_numul(w, 1 / M_norm(w, 2));这一行就指向了一个新生成的matrix,反而是修改后的版本Matrix_copy生成的matrix没有释放,会导致内存泄漏

645770225 commented 3 weeks ago

作者你好,这部分修改代码有一个BUG需要修正一下。 else if (fabs(_x->data[0]) <= 1e-5f) { w = _x; w = M_numul(w, 1 / M_norm(w, 2)); } 中的w=_x应该改为copy,不然_x会在这个函数中提前被free。 应该改成: else if (fabs(_x->data[0]) <= 1e-5f) { w = Matrix_copy(_x); w = M_numul(w, 1 / M_norm(w, 2)); }

好吧,我发现是我本地修改了Mnumul的实现,返回值是新申请了一个矩阵==

645770225 commented 3 weeks ago

不好意思,因为我本地修改了M_numul的实现,返回值是新申请了一个矩阵,所以上面的改动有错误。

  if (_x->data[0] > 1e-5f)
    {
        w = M_add_sub(1, _x, -1, y);
        Matrix *temp_w = w;
        w = M_numul(w, 1 / M_norm(w, 2));
        M_free(temp_w);
    }
    else if (fabs(_x->data[0]) <= 1e-5f)
    {
        w = _x;
        w = M_numul(w, 1 / M_norm(w, 2));
    }
    else
    {
        w = M_add_sub(1, _x, 1, y);
        Matrix *temp_w = w;
        w = M_numul(w, 1 / M_norm(w, 2));
        M_free(temp_w);
    }

应该为

  if (_x->data[0] > 1e-5f)
    {
        w = M_add_sub(1,_x,-1,y);
        M_numul(w,1/M_norm(w, 2));
    }
    else if (fabs(_x->data[0]) <= 1e-5f)
    {
        w = Matrix_copy(_x);
        w = M_numul(w, 1 / M_norm(w, 2));
    }
    else
    {
        w = M_add_sub(1, _x, 1, y);
        M_numul(w, 1 / M_norm(w, 2));
    }