Ask Your Question
1

Kronecker Power of sparse matrices problem

asked 2020-05-02 15:36:12 -0500

SYLA gravatar image

updated 2020-05-12 06:10:53 -0500

Why the following code does not work? I try it also in loops.

import numpy

from scipy import sparse

from scipy.sparse import coo_matrix

def SparseMatrixPower(A,p):
    if p == 1: 
        return(A)
    elif Mod(p,2):
        return(SparseMatrixPower(A,SparseMatrixProduct(A,A)), (p - 1) / 2)  
    else:
        return(SparseMatrixPower(SparseMatrixProduct(A,A), p / 2))

def SparseMatrixProduct(A,B):
        return(sparse.kron(A,B)+sparse.kronsum(A,B))

A=sparse.coo_matrix([[0,1,2],[1,1,2],[2,2,3]])
B=sparse.coo_matrix([[0,1,2],[1,1,2],[2,2,3]])
SparseMatrixProduct(A,B)
SparseMatrixPower(A,3)
edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
1

answered 2020-05-02 17:06:12 -0500

tmonteil gravatar image

updated 2020-05-03 07:43:57 -0500

It seems you are confusing calls to SparseMatrixProduct and SparseMatrixPower. Indeed, SparseMatrixPower(A,A) is meaningless (the second argument should be an integer not a matrix).

When p is odd, your code try to make two steps at once (see the (p-1)/2). To understand the algorithm that is behind your code, I would suggest to separate the add/multiply (p -> p-1) part and the double/squaring (p -> p/2) part as follows :

sage: def SparseMatrixPower(A, p):
....:     if p == 1: 
....:         return(A)
....:     elif Mod(p, 2):
....:         return SparseMatrixProduct(A, SparseMatrixPower(A, p-1))
....:     else:
....:         return SparseMatrixPower(SparseMatrixProduct(A, A), p/2)

That said, if you plan to play a lot with twose objects, I would suggest to define a class for them and define a __mul__method to define the product. Then Sage has generic double-and-add (or multiply-and-square) to automatically define the power, see for example https://doc.sagemath.org/html/en/refe...

edit flag offensive delete link more

Comments

Yes, this is typo, I repaired my code. It does not work as well.

SYLA gravatar imageSYLA ( 2020-05-12 06:10:38 -0500 )edit

It is not completely fixed:

SparseMatrixPower(A,SparseMatrixProduct(A,A))

does not makes sense since the second argument of SparseMatrixPower should be a number, not a matrix.

tmonteil gravatar imagetmonteil ( 2020-05-12 08:31:56 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2020-05-02 15:36:12 -0500

Seen: 86 times

Last updated: May 12