Julia 和 Python 中不同的伪逆矩阵

Question

I am trying to transcribe a python code in Julia. I have a matrix

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

I am computing the (Moore-Penrose) pseudo-inverse of a matrix in Python using numpy.linalg.pinv and the result is

[[-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]]

while in Julia the result of LinearAlgebra.pinv(test) is

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

I would like to ask if anybody knows why the results are different in the two cases and what I can do to have them match. So far I have tried LinearAlgebra.pinv(test[1:3,1:3]), and the result is also different for unknown reasons, but still doesn't match the Python output.

The above 'test' matrix was indeed a test case to simplify the code in a minimum working example, the actual code can be found below.

The full code in Python is:

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)

and the code in Julia is:

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 方面的矩阵构造错误。请注意，range(3)在 python 中与在 Julia中不对应1:3，但是0:2当 python 开始计数时，从零而不是一开始。

更新2：

既然你已经更新了 OP，让我扩展我的答案。如上所述，range(x)在 python 中对应于0:x-1在 Julia 中。同时，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在蟒蛇）。我们对这种差异通过减去1从两个补偿i和j表达式的赋值的RHS（你的具体情况这并不重要，因为移相互补偿）。然后，G_matrix在 Python 和 Julia 中是相同的，并且产生（大致）相同的pinv.