ijv sparse matrix representation for mpmath matrix (coo_matrix equivalent)

Question

Apparently coo_matrix does not work properly for mapmath matrix. Is there any equivalent implementation or simple algorithm to implement for mpmath?

import mpmath
import numpy as np 
from scipy.sparse import coo_matrix

A = np.diag([1]*5)
A_coo = coo_matrix(A)
print(A_coo.row)
print(A_coo.col)
print(A_coo.data)


A_mp = mpmath.diag([1]*5)
A_mp_coo = coo_matrix(A_mp)
print(A_mp_coo.row)
print(A_mp_coo.col)
print(A_mp_coo.data)

output

[0 1 2 3 4]
[0 1 2 3 4]
[1 1 1 1 1]
# ----------------------
[0 0 0 0 0]
[ 0  6 12 18 24]
[mpf('1.0') mpf('1.0') mpf('1.0') mpf('1.0') mpf('1.0')]

Edit:

What I think is loop over every element and find the i,j index of nonzero elements of mpmath matrix then use the following :

size=np.count_nonzero(A_mp)
rows, cols = ... # looping over evey elements to find the indices on nonzero elements

compactMatrix = np.zeros((3, size))
for i in range(size):
    compactMatrix[0, i] = rows[i]
    compactMatrix[1, i] = cols[i]
    compactMatrix[2, i] = A_mp[rows[i], cols[i]]

but it seems vary naive and not fast.

Answer 1

This is what I came up with: I tested it for a square matrix.

first some helper functions:

def to_single_index(i, j, N):
    return N * i + j

def from_single_index(indices, N):
    i_list = []
    j_list = []
    for k in indices:
        i_list.append(k // N)
        j_list.append(k % N)

    return np.array(i_list), np.array(j_list)

def to_mp_matrix(a):

    m, n = a.shape
    A = mpmath.matrix(m,n)
    for i in range(m):
        for j in range(n):
            A[i,j] = a[i, j]
    return A

make a sparse directed matrix to test the code:

N = 5
A = nx.to_numpy_array(nx.gnp_random_graph(
    N, 0.5, seed=2, directed=True))

using coo_matrix to check the correctness:

A_coo = coo_matrix(A)
i = A_coo.row
j = A_coo.col
print(i)
print(j)

# [0 0 1 2 3 3 4 4]
# [3 4 4 4 0 1 2 3]

convert to mpmath matrix and get index of nonzero elements:

A_mp = to_mp_matrix(A)
indices = np.nonzero(A_mp)[0]
print(indices)
# [ 3  4  9 14 15 16 22 23]

test the code on mamath matrix:

indices = to_single_index(i, j, N)
print(indices)
# [ 3  4  9 14 15 16 22 23]

i_list, j_list = from_single_index(indices, N)
print(i_list)
print(j_list)
# [0 0 1 2 3 3 4 4]
# [3 4 4 4 0 1 2 3]