Revision history [back]
 | 1 | initial version |
Tensor product of dense by sparse matrix
Solution in brief
In the example, I and X are (stored as) dense matrices
but Z is (stored as) a sparse matrix.
Strangely, taking the tensor product A.tensor_product(B)
of matrices over a cyclotomic field works when A and B
are (sparse, sparse), (sparse, dense), (dense, dense),
but not (dense, sparse).
Strangely too, the output of the global matrix functions
identity_matrix and diagonal_matrix have different
default sparseness, one returning dense matrices and
the other one sparse matrices.
The sparseness default can be overridden using one of
sparse=True or sparse=False.
One can also use the diagonal_matrix method of the
matrix space where I and X live, instead of the
global diagonal_matrix function.
Below a few words on how to explore the problem, followed by concrete code refactoring suggestions.
Exploration
After the definitions in the question.
Inspect I, X, Z:
sage: I, X, Z
(
[1 0 0] [0 0 1] [ 1 0 0]
[0 1 0] [1 0 0] [ 0 w 0]
[0 0 1], [0 1 0], [ 0 0 -w - 1]
)
Where do they live?
sage: print('\n\n'.join(f"{A.parent()}" for A in (I, X, Z)))
Full MatrixSpace of 3 by 3 dense matrices over ...
Full MatrixSpace of 3 by 3 dense matrices over ...
Full MatrixSpace of 3 by 3 sparse matrices over ...
Try tensor product of I and I: works; of Z and Z: works;
of Z and I: works. But...
Try tensor product of I and Z: fails:
sage: I.tensor_product(Z)
Traceback (most recent call last)
...
AttributeError: 'sage.matrix.matrix_generic_sparse.Matrix_generic_sparse'
object has no attribute '_rational_matrix'
Quick fix:
sage: I.tensor_product(I.parent()(Z))
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 w 0| 0 0 0| 0 0 0]
[ 0 0 -w - 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 1 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 -w - 1| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0| 1 0 0]
[ 0 0 0| 0 0 0| 0 w 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
Code review
Taking this into account, and some simplifications, here are tips for a slight refactoring...
Pick an integer, define a cyclotomic field.
sage: N = 3
sage: k.<w> = CyclotomicField(N)
Optionally, define the matrix space we want to work with:
sage: M = MatrixSpace(k, N)
Two ways to define I:
sage: I = identity_matrix(k, N)
sage: I = M.identity_matrix()
Two ways to define X:
sage: f = lambda i, j: (i - j) % N == 1
sage: X = matrix(k, N, N, f)
sage: X = M(f)
One more way, as a cyclic permutation matrix:
sage: p = Permutation([2 .. N, 1])
sage: X = M(p.to_matrix())
Several ways to define the third matrix:
sage: Z = diagonal_matrix(k, N, [w^j for j in range(N)], sparse=False)
sage: Z = I.parent()(diagonal_matrix([w^j for j in range(N)]))
sage: Z = I.parent().diagonal_matrix([w^j for j in range(N)])
sage: Z = M.diagonal_matrix([w^j for j in range(N)])
Now we have made sure all three matrices are dense.
Their tensor products can be computed. Let us print them all out:
sage: for A, B in ((I, X), (X, I), (I, Z), (Z, I), (X, Z), (Z, X)):
....: print(A.tensor_product(B), end='\n\n\n')
[0 0 1|0 0 0|0 0 0]
[1 0 0|0 0 0|0 0 0]
[0 1 0|0 0 0|0 0 0]
[-----+-----+-----]
[0 0 0|0 0 1|0 0 0]
[0 0 0|1 0 0|0 0 0]
[0 0 0|0 1 0|0 0 0]
[-----+-----+-----]
[0 0 0|0 0 0|0 0 1]
[0 0 0|0 0 0|1 0 0]
[0 0 0|0 0 0|0 1 0]
[0 0 0|0 0 0|1 0 0]
[0 0 0|0 0 0|0 1 0]
[0 0 0|0 0 0|0 0 1]
[-----+-----+-----]
[1 0 0|0 0 0|0 0 0]
[0 1 0|0 0 0|0 0 0]
[0 0 1|0 0 0|0 0 0]
[-----+-----+-----]
[0 0 0|1 0 0|0 0 0]
[0 0 0|0 1 0|0 0 0]
[0 0 0|0 0 1|0 0 0]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 w 0| 0 0 0| 0 0 0]
[ 0 0 -w - 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 1 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 -w - 1| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0| 1 0 0]
[ 0 0 0| 0 0 0| 0 w 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 1 0| 0 0 0| 0 0 0]
[ 0 0 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| w 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 w| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0|-w - 1 0 0]
[ 0 0 0| 0 0 0| 0 -w - 1 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[ 0 0 0| 0 0 0| 1 0 0]
[ 0 0 0| 0 0 0| 0 w 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[--------------------+--------------------+--------------------]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 w 0| 0 0 0| 0 0 0]
[ 0 0 -w - 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 1 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 -w - 1| 0 0 0]
[ 0 0 1| 0 0 0| 0 0 0]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 1 0| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 w| 0 0 0]
[ 0 0 0| w 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[ 0 0 0| 0 0 0|-w - 1 0 0]
[ 0 0 0| 0 0 0| 0 -w - 1 0]
Reparing this in Sage
Fixing this is tracked at:
| | 2 | No.2 Revision |
In the example, I and X are (stored as) dense matrices
but Z is (stored as) a sparse matrix.
Strangely, taking the tensor product A.tensor_product(B)
of matrices over a cyclotomic field works when A and B
are (sparse, sparse), (sparse, dense), (dense, dense),
but not (dense, sparse).
Strangely too, the output of the global matrix functions
identity_matrix and diagonal_matrix have different
default sparseness, one returning dense matrices and
the other one sparse matrices.
The sparseness default can be overridden using one of
sparse=True or sparse=False.
One can also use the diagonal_matrix method of the
matrix space where I and X live, instead of the
global diagonal_matrix function.
Below a few words on how to explore the problem, followed by concrete code refactoring suggestions.
After the definitions in the question.
Inspect I, X, Z:
sage: I, X, Z
(
[1 0 0] [0 0 1] [ 1 0 0]
[0 1 0] [1 0 0] [ 0 w 0]
[0 0 1], [0 1 0], [ 0 0 -w - 1]
)
Where do they live?
sage: print('\n\n'.join(f"{A.parent()}" for A in (I, X, Z)))
Full MatrixSpace of 3 by 3 dense matrices over ...
Full MatrixSpace of 3 by 3 dense matrices over ...
Full MatrixSpace of 3 by 3 sparse matrices over ...
Try tensor product of I and I: works; of Z and Z: works;
of Z and I: works. But...
Try tensor product of I and Z: fails:
sage: I.tensor_product(Z)
Traceback (most recent call last)
...
AttributeError: 'sage.matrix.matrix_generic_sparse.Matrix_generic_sparse'
object has no attribute '_rational_matrix'
Quick fix:
sage: I.tensor_product(I.parent()(Z))
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 w 0| 0 0 0| 0 0 0]
[ 0 0 -w - 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 1 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 -w - 1| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0| 1 0 0]
[ 0 0 0| 0 0 0| 0 w 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
Code review
More variations
Taking this into account, and some simplifications, here are tips for a slight refactoring...
Pick an integer, define a cyclotomic field.
sage: N = 3
sage: k.<w> = CyclotomicField(N)
Optionally, define the matrix space we want to work with:
sage: M = MatrixSpace(k, N)
Two ways to define I:
sage: I = identity_matrix(k, N)
sage: I = M.identity_matrix()
Two ways to define X:
sage: f = lambda i, j: (i - j) % N == 1
sage: X = matrix(k, N, N, f)
sage: X = M(f)
One more way, as a cyclic permutation matrix:
sage: p = Permutation([2 .. N, 1])
sage: X = M(p.to_matrix())
Several ways to define the third matrix:
sage: w_powers = [w^j for j in range(N)]
sage: Z = diagonal_matrix(k, N, [w^j for j in range(N)], w_powers, sparse=False)
sage: Z = I.parent()(diagonal_matrix([w^j for j in range(N)]))
I.parent()(diagonal_matrix(w_powers))
sage: Z = I.parent().diagonal_matrix([w^j for j in range(N)])
I.parent().diagonal_matrix(w_powers)
sage: Z = M.diagonal_matrix([w^j for j in range(N)])
M.diagonal_matrix(w_powers)
Now we have made sure all three matrices are dense.
Their tensor products can be computed. Let us print them all out:
sage: for A, B in ((I, X), (X, I), (I, Z), (Z, I), (X, Z), (Z, X)):
....: print(A.tensor_product(B), end='\n\n\n')
[0 0 1|0 0 0|0 0 0]
[1 0 0|0 0 0|0 0 0]
[0 1 0|0 0 0|0 0 0]
[-----+-----+-----]
[0 0 0|0 0 1|0 0 0]
[0 0 0|1 0 0|0 0 0]
[0 0 0|0 1 0|0 0 0]
[-----+-----+-----]
[0 0 0|0 0 0|0 0 1]
[0 0 0|0 0 0|1 0 0]
[0 0 0|0 0 0|0 1 0]
[0 0 0|0 0 0|1 0 0]
[0 0 0|0 0 0|0 1 0]
[0 0 0|0 0 0|0 0 1]
[-----+-----+-----]
[1 0 0|0 0 0|0 0 0]
[0 1 0|0 0 0|0 0 0]
[0 0 1|0 0 0|0 0 0]
[-----+-----+-----]
[0 0 0|1 0 0|0 0 0]
[0 0 0|0 1 0|0 0 0]
[0 0 0|0 0 1|0 0 0]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 w 0| 0 0 0| 0 0 0]
[ 0 0 -w - 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 1 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 -w - 1| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0| 1 0 0]
[ 0 0 0| 0 0 0| 0 w 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 1 0| 0 0 0| 0 0 0]
[ 0 0 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| w 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 w| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0|-w - 1 0 0]
[ 0 0 0| 0 0 0| 0 -w - 1 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[ 0 0 0| 0 0 0| 1 0 0]
[ 0 0 0| 0 0 0| 0 w 0]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[--------------------+--------------------+--------------------]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 w 0| 0 0 0| 0 0 0]
[ 0 0 -w - 1| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 1 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[ 0 0 0| 0 0 -w - 1| 0 0 0]
[ 0 0 1| 0 0 0| 0 0 0]
[ 1 0 0| 0 0 0| 0 0 0]
[ 0 1 0| 0 0 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 w| 0 0 0]
[ 0 0 0| w 0 0| 0 0 0]
[ 0 0 0| 0 w 0| 0 0 0]
[--------------------+--------------------+--------------------]
[ 0 0 0| 0 0 0| 0 0 -w - 1]
[ 0 0 0| 0 0 0|-w - 1 0 0]
[ 0 0 0| 0 0 0| 0 -w - 1 0]
Reparing this
Summing up
Suggested neat compact version of the code:
sage: N = 3
sage: k.<w> = CyclotomicField(N)
sage: M = MatrixSpace(k, N)
sage: I = M.identity_matrix()
sage: X = M(lambda i, j: (i - j) % N == 1)
sage: Z = M.diagonal_matrix([w^j for j Improving Sage in light of this
Fixing this the bug revealed by the question is tracked at:
There is also a ticket to homogenize the types
of matrix returned by the global functions
identity_matrix and diagonal_matrix
