Conversion between R and Python(indexing issue?)-CodePudding

I'm trying to convert this R function to Python. However, I am running into issues with the results being incorrect.

The R function in question:

construct_omega <- function(k){

    E <- diag(2*k)
    omega <- matrix(0, ncol=2*k, nrow=2*k)

    for (i in 1:k){
        omega <- omega  
            E[,2*i-1] %*% t(E[,2*i]) -
            E[,2*i] %*% t(E[,2*i-1])
    }

    return(omega)

}

This is my current attempt at porting the function to Python:

def construct_omega(k=1):

    E = np.identity(2*k)
    omega = np.zeros((2*k, 2*k))

    for i in range(1, k):
        omega = omega   \
            E[:,2*i-1] * np.transpose(E[:,2*i]) - \
            E[:,2*i] * np.transpose(E[:,2*i-1])
    
    return omega

In R, the result matrix is this:

> construct_omega(2)
     [,1] [,2] [,3] [,4]
[1,]    0    1    0    0
[2,]   -1    0    0    0
[3,]    0    0    0    1
[4,]    0    0   -1    0

But in Python, the result is the 4x4 zero matrix.

Any help would be appreciated, thanks!

CodePudding user response：

The issue here is an edge case of matrix multiplication in numpy. You should consult the docs and this post or this one. Basically what happens is that you are doing a dot product and getting a scalar, not a matrix as you assume. The fix is to use np.outer. The other issue is the indexing that in python starts at 0 so you need to rewrite your code a bit.

import numpy as np

def construct_omega(k=1):
    
    E = np.identity(2*k)
    omega = np.zeros((2*k, 2*k))

    for i in range(k):
        omega = omega   \
            np.outer(E[:,2*i], E[:,2*i 1].T) - \
            np.outer(E[:,2*i 1], E[:,2*i].T)
        
    return omega

construct_omega(1)
#array([[ 0., 1.],
#       [ -1.,  0.]])

and

construct_omega(2)
#array([[ 0.,  1.,  0.,  0.],
#       [-1.,  0.,  0.,  0.],
#       [ 0.,  0.,  0.,  1.],
#       [ 0.,  0., -1.,  0.]])