Julia中相同矩阵乘积结果不一致的问题排查与修正求助

在Jupyter Lab运行Julia代码时，执行相同的矩阵乘积操作，两次循环输出的final_vec结果差异明显。代码如下：

using LinearAlgebra

function U1(α, M)
    return exp(α * M)
end

function U2(β, M)
    return exp(β * M)
end

test_list=[1 2]

α_list = test_list[1]
β_list = test_list[2]

M1=[1 3 3; 3 2 2; 2 2 1]
M2=[3 2 1; 3 0 3; 2 1 1]

init_vec=[1 ; 1 ; 1]

for _ in 1:2
    α=α_list[1]
    β=β_list[1]
    final_vec = U1(α,M1) * U2(β,M2) * init_vec
    println(final_vec)
end

# 输出结果：
# [3.2498540882389616e7, 3.362948386386128e7, 2.2549887512082998e7]
# [3.3668507705110244e7, 3.4777346320714206e7, 2.548964322110533e7]

问题成因分析 #

• 交互式环境变量污染：Jupyter Lab的交互式特性可能导致之前运行的单元格意外修改了M1、M2、α_list或β_list的值，使得两次循环使用的参数不一致。
• 索引逻辑潜在错误：代码中α_list = test_list[1]和β_list = test_list[2]的索引方式存在歧义：若test_list是1×2矩阵，test_list[1]会取到标量1，后续α_list[1]会触发语法错误；若实际运行时test_list是动态数组（如未固定的随机结构），可能导致每次循环取到的α/β值不同。
• 极小概率的数值计算非确定性：使用多线程BLAS库时，矩阵运算的顺序差异可能导致极小数值偏差，但你遇到的结果差异过大，这种可能性极低。

修改方案 #

方案1：简化代码并固定变量定义

去掉多余的α_list/β_list，直接使用原始变量，同时确保所有依赖变量在单元格内重新定义，避免环境污染：

using LinearAlgebra

function U(θ, M)
    return exp(θ * M)
end

# 固定所有输入变量，放在单元格最顶部
test_list = [1, 2]  # 使用一维向量而非矩阵，避免索引歧义
M1 = [1 3 3; 3 2 2; 2 2 1]
M2 = [3 2 1; 3 0 3; 2 1 1]
init_vec = [1; 1; 1]

for _ in 1:2
    α = test_list[1]
    β = test_list[2]
    # 打印参数确认一致性
    println("当前参数：α=$α, β=$β")
    final_vec = U(α, M1) * U(β, M2) * init_vec
    println(final_vec)
end

方案2：验证参数一致性

在循环内添加参数打印，确认每次循环使用的矩阵和标量完全一致：

for _ in 1:2
    α = α_list[1]
    β = β_list[1]
    println("α: $α, β: $β")
    println("M1:\n$M1")
    println("M2:\n$M2")
    final_vec = U1(α,M1) * U2(β,M2) * init_vec
    println(final_vec)
end

如果发现参数不一致，排查是否有其他单元格修改了这些变量。

方案3：使用确定性计算

若怀疑矩阵指数运算的稳定性，可手动实现小矩阵的泰勒展开（仅适用于小矩阵）：

function mat_exp(M, n_terms=20)
    result = Matrix{Float64}(I, size(M)...)
    current_term = result
    for k in 1:n_terms
        current_term = current_term * M / k
        result += current_term
    end
    return result
end

# 修改U1/U2使用自定义的mat_exp
function U1(α, M)
    return mat_exp(α * M)
end

function U2(β, M)
    return mat_exp(β * M)
end

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