如何用Gonum库实现对标Matlab chol()的Cholesky分解？

问题分析 #

你之前的代码存在两个核心问题：

• 输入的对称矩阵与Matlab中的P_inno完全不符，导致分解结果毫无意义
• 对UTo方法的用法理解有误，且未正确输出分解后的矩阵

修正后的完整代码 #

package main

import (
	"fmt"
	"gonum.org/v1/gonum/mat"
)

func main() {
	// 按gonum SymDense的要求构造Matlab中的P_inno矩阵
	// SymDense的data是下三角元素按列主序排列：
	// 第0列: (0,0), (1,0), (2,0), (3,0)
	// 第1列: (1,1), (2,1), (3,1)
	// 第2列: (2,2), (3,2)
	// 第3列: (3,3)
	arr := []float64{
		0.10062,
		-0.042635, 0.21293,
		-0.072741, -0.02717, 0.26536,
		-0.1434, -0.052181, 0.27184, 0.86335,
	}
	P_inno := mat.NewSymDense(4, arr)

	// 执行Cholesky分解
	var chol mat.Cholesky
	if ok := chol.Factorize(P_inno); !ok {
		fmt.Println("矩阵不是正定/半正定矩阵")
		return
	}

	// 创建TriDense实例存储上三角分解结果
	var t mat.TriDense
	// 将分解结果写入t，检查操作是否成功
	if err := chol.UTo(&t); err != nil {
		fmt.Println("提取分解结果失败:", err)
		return
	}

	// 格式化输出，与Matlab风格对齐
	fmt.Println("chol(P_inno) = ")
	fmt.Printf("%v\n", mat.Formatted(&t, mat.Prefix(""), mat.Squeeze()))
}

关键要点说明 #

• 对称矩阵初始化：SymDense要求传入的data数组是矩阵的下三角元素按列主序排列，不能用你之前的对角矩阵填充，这是匹配Matlab输入的核心前提。
• UTo方法的正确使用：
  • UTo的作用是将Cholesky分解得到的上三角矩阵写入传入的*TriDense变量，它返回的是error（用于检查操作异常），而非矩阵本身。
  • 必须传入&t（指针），而不是t（值传递），否则无法修改t的内容。
• 矩阵输出：直接打印TriDense变量无法得到友好的矩阵格式，需要用mat.Formatted进行格式化，通过mat.Prefix("")和mat.Squeeze()参数让输出与Matlab的风格一致。
• 分解逻辑：chol.Factorize(P_inno)会修改chol变量本身，将分解结果存储在其中，返回的布尔值仅用于判断矩阵是否可分解（正定/半正定），这是符合gonum库设计的正常逻辑。

输出结果 #

运行上述代码后，输出将与Matlab的chol(P_inno)结果一致：

chol(P_inno) = 
0.31721  -0.13441  -0.22932  -0.45207
0        0.44143  -0.13137  -0.25585
0        0        0.44217   0.30432
0        0        0         0.70775

内容的提问来源于stack exchange，提问作者Chrissy Heu