Search code examples
pythonmatlabscipysparse-matrix

Translate matlab code to python (scipy)


I' trying to translate this matlab code to python:

T = length(z);
lambda = 10;
I = speye(T)
D2 = spdiags(ones(T-2,1)*[1 -2 1],[0:2],T-2,T);
z_stat = (I-inv(I + lambda^2*D2'*D2))*z;

What I got at the moment:

T = len(signal)
lam = 10;
I  = np.identity(T)
D2 = scipy.sparse.spdiags(np.ones((T-2,1),dtype=np.int)*[1,-2,1],(range(0,3)),T-2,T);

At the moment I get this error

"scipy.sparse.sp...ge(0,3)),T-2,T)" ValueError: number of diagonals (298) does not match the number of offsets (3) args tuple: ('number of diagonals (298) does not match the number of offsets (3)',)

When looking at the documentation, the matlab function and the python function are very similar. Though there is probably a difference which I am missing. My question is now: What am I doing wrong ?

edit: z is an array with length 300


Solution

  • If you transpose the filter data to spdiags then you get an answer of the same dimensions in both packages:

    # numpy/scipy
    filt = [1,-2,1]* np.ones((1,T-2),dtype=np.int).T
    D2 = scipy.sparse.spdiags(data.T, (range(0,3)),T-2,T)
    np.shape(D2)
    >>> (298, 300)
    
    % matlab check
    D2 = spdiags(ones(T-2,1)*[1 -2 1],[0:2],T-2,T)
    size(D2)
    ans =
      298   300