Julia和Python中的伪逆矩阵不同

Question

我正在尝试在 Julia 中转录 Python 代码。我有一个矩阵

test = [2.0 3.0 4.0
        3.0 4.0 5.0
        4.0 5.0 6.0]

我正在使用numpy.linalg.pinv 在 Python 中计算矩阵的 (Moore-Penrose) 伪逆，结果是

[[-8.33333333e-01 -1.66666667e-01  5.00000000e-01]
 [-1.66666667e-01 -7.86535165e-17  1.66666667e-01]
 [ 5.00000000e-01  1.66666667e-01 -1.66666667e-01]]

在 Julia 中，LinearAlgebra.pinv(test) 的结果是

3×3 Array{Float64,2}:
 -1.33333   -0.166667      1.0     
 -0.166667  -6.59195e-17   0.166667
  1.0        0.166667     -0.666667

我想问是否有人知道为什么这两种情况下的结果不同，以及我可以做些什么来使它们匹配。到目前为止，我已经尝试过LinearAlgebra.pinv(test[1:3,1:3])，结果也因未知原因而有所不同，但仍然与Python输出不匹配。

上面的“测试”矩阵确实是一个测试用例，以简化最小工作示例中的代码，实际代码可以在下面找到。

Python中的完整代码是：

import numpy as np
import random as rd
import matplotlib.pyplot as plt

n_bin = 8
A = 5.
sigma_G = 3.

G_temp = np.zeros((n_bin,n_bin))

for i in range(n_bin):
     for j in range(n_bin):
         G_temp[i,j] = A*np.exp(-(1/2)*((i-j)**2)/sigma_G**2)

G_matrix = np.matmul(G_temp,G_temp.T)
G_inv = np.linalg.pinv(G_matrix)

Julia 中的代码是：

using LinearAlgebra
using Distributions  
using Base

n_bin = 8
A = 5. 
sigma_G = 3. 

G_temp = zeros(n_bin,n_bin)

for i = 1:n_bin
    for j = 1:n_bin
         G_temp[i,j] = A*exp(-(1/2)*((i-j)^2)/sigma_G^2)
    end
end

G_matrix = G_temp*transpose(G_temp)
G_inv = LinearAlgebra.pinv(G_matrix)

Answer 1

我无法重现：

julia> using PyCall; np = pyimport("numpy");

julia> test = [2. 3. 4.; 3. 4. 5.; 4. 5. 6.]
3×3 Array{Float64,2}:
 2.0  3.0  4.0
 3.0  4.0  5.0
 4.0  5.0  6.0

julia> pinv(test)
3×3 Array{Float64,2}:
 -1.33333   -0.166667      1.0
 -0.166667  -4.16334e-17   0.166667
  1.0        0.166667     -0.666667

julia> using PyCall; np = pyimport("numpy");

julia> pinv(test) ≈ np.linalg.pinv(test)
true

请注意，与您发布的内容相比，我使用 numpy 得到了不同的伪逆。

>>> test = [[2.,3.,4.],[3.,4.,5.],[4.,5.,6.]]
>>> import numpy as np
>>> np.linalg.pinv(test)
array([[ -1.33333333e+00,  -1.66666667e-01,   1.00000000e+00],
       [ -1.66666667e-01,  -2.42861287e-17,   1.66666667e-01],
       [  1.00000000e+00,   1.66666667e-01,  -6.66666667e-01]])

更新：

您发布的 numpy 输出对应于以下略有不同的矩阵的伪逆：

julia> test2 = [0. 1. 2.; 1. 2. 3.; 2. 3. 4.]
3×3 Array{Float64,2}:
 0.0  1.0  2.0
 1.0  2.0  3.0
 2.0  3.0  4.0

julia> pinv(test2)
3×3 Array{Float64,2}:
 -0.833333  -0.166667      0.5
 -0.166667  -7.63278e-17   0.166667
  0.5        0.166667     -0.166667

这可能是 python 端的矩阵构造错误。请注意，python 中的range(3) 确实不 对应于 Julia 中的 1:3，但 0:2 因为 python 从零而不是一开始计数。

更新 2：

既然你已经更新了 OP，让我扩展我的答案。如上所述，python 中的range(x) 对应于 Julia 中的0:x-1。同时，python 中的数组索引从 0 开始，而在 Julia 中它从 1 开始。因此，假设 python 代码产生“正确”（预期）结果，您的 Julia 代码应如下所示：

using LinearAlgebra
# dropped Base and Distributions here since you don't need them.

n_bin = 8
A = 5. 
sigma_G = 3. 

G_temp = zeros(n_bin,n_bin)

for i = 1:n_bin
    for j = 1:n_bin
         G_temp[i,j] = A*exp(-(1/2)*(((i-1)-(j-1))^2)/sigma_G^2) # note the (i-1) and (j-1) here!
    end
end

G_matrix = G_temp*transpose(G_temp)
G_inv = pinv(G_matrix)

请注意循环变量i 和j，用于索引G_temp 从1 到n_bin（而不是0 到n_bin-1，如在python 中）。我们通过在右上方的表达式中从 i 和 j 中减去 1 来补偿这种差异。分配（在您的特定情况下，这无关紧要，因为班次相互补偿）。然后，G_matrix 在 Python 和 Julia 中是相同的，并且产生（大致）相同的 pinv。