Why are coefficients in tensor product bases represented the way they are?

If you have a basis of Chebyshev Polynomials, for example, that looks like this:

[T00 T01 T02;
 T10 T11  0 ;
 T20  0   0 ]

These are stored like so

[T00, T01, T10, T02, T11, T20]

My question is: why? Would it not be more natural to simply keep using the matrix? Or perhaps reshape the matrix to a vector, so columns are concatenated together?

In my case, using the matrix directly means I don't have a differentiation operator with 1e12 entries, but rather one with 1e6 entries, which is the difference between feasible and impossible. (The reduction occurs because in the matrix way of doing things, I can use the 1D operators on each column/row, which amounts to a simple matrix matrix multiplication).

Of course, this comes at the cost of doubling the memory, but it seems more than worth it to double the memory at this stage if you can save many orders of magnitude in memory/time when building operators and applying them.

No files, tests, or entry points are named. Start by locating the tensor-product basis coefficient representation and the differentiation-operator construction, then compare the current packed ordering with matrix or reshaped storage for the Chebyshev example. Done means documenting the rationale or agreeing on a feasible representation change that addresses the reported memory and performance costs.