Zero-dimensional matrices

Question

In a recent Mathematica project, $(n \times 0)$- and $(0 \times n)$-dimensional matrices have suddenly become a frustratingly common edge case for me. For instance, consider the following two situations, both common in linear algebra programming:

1. Suppose we have a matrix $A$, presenting a linear transformation $\newcommand{\R}{\mathbb{R}}\R^n \to \R^m$. The kernel of $A$ is a particular subspace of the source, $\R^n$, and if $A$ has a kernel of dimension $d > 0$, then NullSpace[A] produces a new matrix presenting a map $\R^d \to \R^n$ whose image is $\ker A$. This mode of thought is very useful; it lets you, for instance, restrict the action of some other map $B: \R^n \to \R^l$ to the kernel of $A$ by multiplication: B.Transpose[NullSpace[A]].

2. Suppose we have a pair of matrices $A$ and $B$, presenting a pair of linear transformations $\R^n \xrightarrow{A} \R^l \xleftarrow{B} \R^m$ with a common target. The pullback or generalized intersection of these two maps is defined to be the collection of pairs of vectors $(x, y) \in \R^n \times \R^m$ with $Ax = By$. This, too, is easy enough to model in Mathematica: the NullSpace of the block matrix $\left[\begin{array}{c|c}A & -B\end{array}\right]$ splits into a block matrix itself $\left[\begin{array}{c}P \\ \hline Q \end{array}\right]$. The matrices $P$ and $Q$ present the "largest nontrivial" maps into $\R^n$ and $\R^m$ respectively with the property that $AP = BQ$.

While the mathematics continues to make sense, both of these run into serious difficulties when $0$-dimensional vector spaces get involved. This happens all the time in both situations (and, generally, in any situation where you might be applying NullSpace to a matrix with trivial kernel): because Mathematica models matrices as lists of lists (with some exceptions, like SparseArrays, which have to support conversion to lists of lists anyway), it has real trouble modeling an $(n \times 0)$-dimensional matrix. For instance, NullSpace[ IdentityMatrix[n]] returns {}, which at first glance is reasonable enough but is totally unsuited for uniformly handling the above procedures. The expression {} is not $(n \times 0)$-dimensional in any sense, and so matrix multiplication, transposition(, inversion in the case $n = 0$), and so on all balk at using this as input.

This is surmountable, of course, by riddling my Mathematica code with conditionals that check for this failure output and handle it separately. This is uglier than it sounds (and it already sounds ugly); oftentimes this means passing around extra information to keep track of the n involved, since it's not recoverable from {}.

My question is: Is there a cleaner solution than a sea of conditionals? Is there an idiomatic way to model $(n \times 0)$- and $(0 \times n)$-dimensional matrices in Mathematica, ideally in a way that's near-seamlessly compatible with standard routines like Times, Transpose, ArrayFlatten, ...?

(P.S.: I eventually intend to offer this thing I'm writing as a package, and so I'd prefer not to modify the behavior of built-in functions directly.)

Answer 1

A little late, I guess, since I had similar ideas to both rm-rf and ybeltukov.

---

I would try using one head to represent an empty matrix, and define what Dot etc. means via TagSetDelayed. You can add error messages for mismatched dimensions in the Condition test -- you'll certainly need them.

ClearAll[emptyMatrix];

  (* this comes up in the product of an r x 0 times 0 x c empty matrices *)
emptyMatrix[r_Integer, c_Integer] /; r > 0 && c > 0 := ConstantArray[0, {r, c}];

 (* extend Dot *)
emptyMatrix /: 
  Dot[m1_emptyMatrix, m2_?MatrixQ] /; 
   Last[m1] == First@Dimensions[m2] := 
  emptyMatrix[First@Dimensions[m1], Last[m2]];
emptyMatrix /: 
  Dot[m1_?MatrixQ, m2_emptyMatrix] /; 
   Last@Dimensions[m1] == First[m2] := 
  emptyMatrix[First@Dimensions[m1], Last[m2]];
emptyMatrix /: 
  Dot[m1_emptyMatrix, m2_emptyMatrix] /; Last[m1] == First[m2] := 
  emptyMatrix[First[m1], Last[m2]];

  (* extend Transpose *)
emptyMatrix /: Transpose[emptyMatrix[r_, c_]] := emptyMatrix[c, r];

  (* extend scalar multiplication -- other uses of Times left as an exercise *)
emptyMatrix /: Times[k_?NumericQ, m_emptyMatrix] := m;

(* etc. *)

I think you have to define your own interfaces to some linear algebra functions. For example, NullSpace:

Clear[nullSpace];
nullSpace[m_] := NullSpace[m] /. {} -> emptyMatrix[First@Dimensions[m], 0];

ArrayFlatten seems like it might take some work to get it right.

Examples:

nullSpace[IdentityMatrix[3]]
(*
   emptyMatrix[3, 0]
*)

RandomReal[1, {3, 3}].nullSpace[IdentityMatrix[3]]
(*
   emptyMatrix[3, 0]
*)

RandomReal[1, {2, 3}].nullSpace[IdentityMatrix[3]]
(*
   emptyMatrix[2, 0]
*)

emptyMatrix[4, 0].emptyMatrix[0, 3]
(*
   {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}
*)

{{}}.emptyMatrix[0, 3]
(*
   {{0, 0, 0}}
*)

{{}, {}, {}}.emptyMatrix[0, 3]
(*
   {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}
*)

RandomReal[1, {3, 2}].nullSpace[IdentityMatrix[3]]
(* Here an error message would be nice
   {{0.0973613, 0.422063}, {0.896662, 0.522105}, {0.357575, 0.980591}}.emptyMatrix[3, 0]
*)

Transpose[emptyMatrix[0, 3]]
(*
   emptyMatrix[3, 0]
*)

2 emptyMatrix[5, 0]
(*
   emptyMatrix[5, 0]
*)

Other problems might include extending Map and Fold.

---

I can add to rm-rf's answer that ArrayFlatten works with $n \times 0$ matrices of the form {{}, {}, ..}, but you're still left with the problem of $0 \times n$ matrices.

m1 = {{1, 2, 3}, {3, 4, 5}};
m2 = {{1, 2}, {3, 4}};
m0 = {{}, {}};

ArrayFlatten[{{0, 0, m1}, {m2, m0, 0}}] // MatrixForm

This might be helpful in working out how to extend ArrayFlatten to general empty matrices.

Answer 2

Another possibility is definition of your own types

ClearAll[nxo, oxn];
nxo /: Transpose[nxo[n_]] := oxn[n];
oxn /: Transpose[oxn[n_]] := nxo[n];
nxo /: A_ .nxo[n_] /; Dimensions[A][[2]] == n := nxo@Dimensions[A][[1]];
oxn /: oxn[n_].A_ /; Dimensions[A][[1]] == n := oxn@Dimensions[A][[2]];

Note the space between A_ and .nxo. You can add much more rules for another built-in commands.

IdentityMatrix[{5, 10}].nxo[10]

nxo[5]

Answer 3

You can create $n\times 0$ matrices in Mathematica as follows:

nx0[n_Integer] := ConstantArray[{}, n]
Dimensions@nx0[10]
(* {10, 0} *)

This will also obey matrix multiplication (and perhaps a few other operations):

IdentityMatrix[10].nx0[10]
(* {{}, {}, {}, {}, {}, {}, {}, {}, {}, {}} *)

However, I don't know of a good way to natively represent $0\times n$ matrices. Transposing the above will collapse all the nested empty lists to a single {}, which is not useful.

This doesn't go all the way for what you wanted, but perhaps it can be of some use.