Pseudo-inverse matrix different in Julia and 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

I can not reproduce:

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

Note that I get a different pseudo inverse using numpy compared to what you posted.

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

UPDATE:

The numpy output you posted corresponds to the pseudo inverse of the following slightly different matrix:

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

This is likely a matrix construction mistake on the python side. Note that range(3) in python does not correspond to 1:3 in Julia, but 0:2 as python starts counting with zero rather than one.

UPDATE2:

Since you have updated the OP, let me extend my answer. As pointed out above, range(x) in python corresponds to 0:x-1 in Julia. At the same time, array indexing in python starts at 0 while in Julia it starts with 1. Hence, assuming the python code produces the "correct" (expected) result, your Julia code should be as follows:

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)

Note that the loop variables i and j, which are used for indexing G_temp go from 1 to n_bin (and not 0 to n_bin-1 as in python). We compensate for this difference by substracting 1 from both i and j in the expression on the r.h.s. of the assignment (in your particular case this doesn't matter because the shifts compensate each other). Then, the G_matrix is the same in Python and Julia and produces (roughly) the same pinv.