M_eigen求特征向量有误

作者您好，感谢开源测试发现，M_eigen求特征向量时偶尔出现不对的情况，如：

    float test_data[9]={3.0f,2.0f,3.0f,4.0f,5.0f,6.0f,2.0f,8.0f,9.0f};
    Matrix *test=Matrix_gen(3,3,test_data);
    Matrix **eig = M_eigen(test);
    printf("orign matrix\n");
    M_print(test);
    printf("eigen_vel\n");
    M_print(eig[0]);
    printf("eigen_vec\n");
    M_print(eig[1]);
    printf("orign*eigen_vec\n");
    M_print(M_mul(test,eig[1]));

    Matrix * vel_diag=M_I(test->row);
    for(int i =0;i<vel_diag->row;i++)
    {
        vel_diag->data[i*vel_diag->column+i]=eig[0]->data[i];
    }
    printf("eigen_vec\n");
    M_print(vel_diag);
    printf("vel * eigen_vec\n");
    M_print(M_mul(eig[1],vel_diag));

计算结果为： orign matrix 3.000 2.000 3.000
4.000 5.000 6.000
2.000 8.000 9.000
eigen_vel 15.416 1.000 0.584
eigen_vec 0.279 -0.333 -0.323
0.557 -0.667 -0.646
0.782 0.667 0.691
orign*eigen_vec 4.297 -0.333 -0.189
8.594 -0.667 -0.377
12.055 -0.000 0.404
eigen_vec 15.416 0.000 0.000
0.000 1.000 0.000
0.000 0.000 0.584
vel * eigen_vec 4.297 -0.333 -0.189
8.594 -0.667 -0.377
12.055 0.667 0.404

算出的原矩阵*特征向量与特征向量*对角化的特征值矩阵结果不一，特征值为1时的特征向量计算不对，检查发现是由特征值求特征向量过程中，（A-lamda*I）矩阵化为阶梯型矩阵时，当最后一行与前面行线性相关时才能计算正确的特征向量，此处原理实质上是高斯消元法，却没有判断主元以及主元所在行与非主元行的交换操作，因为（A-lamda*I）必定存在两行线性相关，当线性相关的行不在最后一行时，就会出现错误，此外，当主元为0时，也会出现计算错误，因此，我将M_eigen改为：

Matrix **M_eigen(Matrix *_mat)
{
    Matrix **M_array_eigen_vec = NULL;
    if (_mat->column == _mat->row)
    {
        M_array_eigen_vec = (Matrix **)malloc(sizeof(Matrix *) * 2); // 保存Q/R矩阵地址
        enum
        {
            val = 0,
            vec = 1
        };
        Matrix *eigen_value = M_eigen_val(_mat);
        float temp;
        M_array_eigen_vec[val] = eigen_value;
        int eigen_count, dim = _mat->column, i, j, k, ik, jk;
        Matrix *eigen_vector = NULL, *_mat_ = NULL;
        eigen_vector = M_Zeros(dim, dim); // 生成特征向量
        M_array_eigen_vec[vec] = eigen_vector;
        MATRIX_TYPE eigen_value_temp;
        MATRIX_TYPE coe; // core of elements, 对角线元素/中心元素
        for (eigen_count = 0; eigen_count < dim; eigen_count++)
        {
            _mat_ = Matrix_copy(_mat);
            eigen_value_temp = eigen_value->data[eigen_count];

            for (i = 0; i < dim; i++)
            {
                _mat_->data[i * _mat_->column + i] -= eigen_value_temp; // 注意: 这里计算 (A-lamda*I), 当A为I/diag时，可能存在问题;
            }

            // 矩阵化为阶梯型矩阵(归一性): 对角线值为1
            for (i = 0; i < dim - 1; i++)
            {
                coe = _mat_->data[i * dim + i];
                k = i;
                for (j = i + 1; j < dim; j++)//寻找主元
                {
                    if (fabs(_mat_->data[j * dim + i]) > fabs(coe))
                    {
                        coe = _mat_->data[j * dim + i];
                        k = j;
                    }
                }
                if (fabs(coe) < 1e-5f) // 如果主元为0
                {
                    continue;
                }
                if (k != i)//交换主元与非主元行
                {
                    for (j = 0; j < dim; j++)
                    {
                        temp = _mat_->data[i * dim + j];
                        _mat_->data[i * dim + j] = _mat_->data[k * dim + j];
                        _mat_->data[k * dim + j] = temp;
                    }
                }

                for (j = i; j < dim; j++)
                {
                    _mat_->data[i * dim + j] /= coe; // 让对角线元素归一化
                }

                for (ik = i + 1; ik < dim; ik++)
                {
                    coe = _mat_->data[ik * dim + i];
                    for (jk = i; jk < dim; jk++)
                    {
                        _mat_->data[ik * dim + jk] -= coe * _mat_->data[i * dim + jk];
                    }
                }
            }
            // printf("eigen_value_temp:%f\n", eigen_value_temp);
            // M_print(_mat_);
            MATRIX_TYPE sum1 = 1;
            if (abs(_mat_->data[(dim - 1) * dim + (dim - 1)]) > 1e-5f) // 如果矩阵最后一行不为0
            {
                sum1 = 0;
                eigen_vector->data[(dim - 1) * dim + eigen_count] = 0.0f;
            }
            else
            {
                sum1 = 1;
                eigen_vector->data[(dim - 1) * dim + eigen_count] = 1;
            }
            for (ik = dim - 2; ik >= 0; ik--)
            {
                MATRIX_TYPE sum2 = 0;
                for (jk = ik + 1; jk < dim; jk++)
                {
                    sum2 += _mat_->data[ik * dim + jk] * eigen_vector->data[jk * dim + eigen_count];
                }
                if (fabs(_mat_->data[ik * dim + ik]) > 1e-5f) // 如果对角线元素不为0
                {
                    sum2 = -sum2 / _mat_->data[ik * dim + ik];
                }
                else
                {
                    sum2 = 1;
                }
                sum1 += sum2 * sum2;
                eigen_vector->data[ik * dim + eigen_count] = sum2;
            }
            M_free(_mat_);
            sum1 = sqrt(sum1); // 当前列的模
            int i;
            for (i = 0; i < dim; i++)
            {
                // 向量单位化
                eigen_vector->data[i * dim + eigen_count] /= sum1;
            }
        }
    }
    else
    {
        return NULL;
    }
    return M_array_eigen_vec;
}

计算结果为： orign matrix 3.000 2.000 3.000
4.000 5.000 6.000
2.000 8.000 9.000
eigen_vel 15.416 1.000 0.584
eigen_vec 0.279 -0.456 -0.323
0.557 -0.570 -0.646
0.782 0.684 0.691
origneigen_vec 4.297 -0.456 -0.189
8.594 -0.570 -0.377
12.055 0.684 0.404
eigen_vec 15.416 0.000 0.000
0.000 1.000 0.000
0.000 0.000 0.584
vel eigen_vec 4.297 -0.456 -0.189
8.594 -0.570 -0.377
12.055 0.684 0.404

可以看到，特征向量计算已经正常

Amoiensis / Matrix_hub

M_eigen求特征向量有误 #10