# 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]
```

## 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_powers, sparse=False)
sage: Z = I.parent()(diagonal_matrix(w_powers))
sage: Z = I.parent().diagonal_matrix(w_powers)
sage: Z = 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]
```

## 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 in range(N)])
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')
```

## Improving Sage in light of this

Fixing 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`