Is there a (matrix) operation that can count elements in vector pairs?

For two vectors $$x \in \{0, 1, 2\}^{n}$$ and $$y \in \{0, 1, 2\}^{n},$$ I need to generate a matrix $$C\in \mathcal{R}^{3\times3}$$ where $$C_{i,j}$$ equals the number of indexes $$t$$ where $$x[t]=i$$ and $$y[t]=j$$.

I know I can always write a program to iteratively check and count, but I wonder if there is any way that I can approximate this with matrix operations?

Or, is there any way that I can approximate this as an optimization problem?

• I don't quite understand what you're asking. Are you looking for mathematical notation for how to express $C$ given $x$ and $y$? Or do you want a way to calculate it? Also, I don't understand why approximation and optimization are mentioned, as they don't really seem to be appropriate/necessary here given what you've said so far. – user20160 Aug 3 '18 at 9:04
• @user20160 I'm looking for how to calculate it (other than iteratively checking and counting). Since there is probably no way to calculate it directly, approximation through some optimization is also fine. – Haohan Wang Aug 3 '18 at 13:41
• Could you explain what you mean by "iteratively check and count"? After all, all you're doing is tabulating the nine possible values of $(x[t], y[t]).$ No iteration or checking is ordinarily needed for such a straightforward and simple calculation. – whuber Aug 3 '18 at 15:00

It comes down to what you consider a "matrix operation." Computer programming platforms usually admit a larger gamut of possibilities than linear algebra, for instance. Let's focus on the former. (There are solutions using the most basic linear algebraic operations, but they are a bit more involved: the indicator function can be implemented as a polynomial.)

One solution is afforded by a generalized outer product. Given any binary function $$f$$ and two vectors $$\mathbf{x}$$ and $$\mathbf{y},$$ the outer product $$\mathbf{x}\, {\otimes}_f\ \mathbf{y}$$ is the matrix

$$\left(\mathbf{x}\, {\otimes}_f\ \mathbf{y}\right)_{ij} = f(x_i, y_j).$$

Let $$s$$ be any scalar (a potential component of a vector). The indicator of $$s,$$ written $$I(,\, s),$$ is the function that returns the value $$1$$ when applied to $$s$$ and otherwise returns $$0:$$

$$I(t,s) = \cases{1 & t=s \\0 & \text{otherwise.}}$$

Suppose the $$m$$-vector $$\mathbf{x}$$ has values in $$\{0,1,\ldots, n\}.$$ The outer product of $$x$$ with the vector $$\mathbf{n} = (0,1,2,\ldots,n)$$ with respect to the indicator function $$I$$ is the $$m\times n$$ matrix

$$\mathbf{x}\, {\otimes}_I\, \mathbf{n}$$

whose $$i,j$$ entry is the indicator $$I(x_i, j-1).$$ When $$\mathbf{x}$$ and $$\mathbf{y}$$ are both $$m$$-vectors, the matrix product

$$\left(\mathbf{x}\, {\otimes}_I\, \mathbf{n}\right)^\prime \left(\mathbf{y}\, {\otimes}_I\, \mathbf{n}\right) = \left(\mathbf{n}\, {\otimes}_I\, \mathbf{x}\right) \left(\mathbf{y}\, {\otimes}_I\, \mathbf{n}\right)$$

Proof For $$0 \le i, j\le n,$$ just calculate from the foregoing definitions that

$$\left(\left(\mathbf{x}\, {\otimes}_I\, \mathbf{n}\right)^\prime \left(\mathbf{y}\, {\otimes}_I\, \mathbf{n}\right)\right)_{i+1,j+1} = \sum_{t=1}^n I(x_t, i) I(y_t, j).$$

The terms in the sum are zero except where $$x_t=i$$ and $$y_t=j,$$ where they are equal to $$(1)(1)=1.$$ The sum therefore counts all such terms, which is exactly what is intended.

The outer products each require $$O(mn)$$ computational effort, while the multiplication of the resulting $$n\times m$$ by $$m\times n$$ matrices takes $$O(n^2m)$$ effort and therefore dominates the cost. This appears inefficient, because the same result can be obtained simply by tabulating the values by (a) initializing the final $$n\times n$$ array and (b) scanning once across both vectors simultaneously, updating the counts in the array as you go. That is only $$O(n^2+m)$$ effort.

The interest in this approach therefore would focus either on being able to apply algebraic rules to analyze the operation or to exploit built-in efficiencies on an array-oriented computing platform. At the end of this post I assess the latter possibility for R.

Here is an implementation in R.

f <- function(x, y, n=2) outer(0:n, x, ==) %*% outer(y, 0:n, ==)

As an example of its use, let's generate a pair of vectors exhibiting many possible combinations of values from $$\{0,1,2\}$$:

X <- expand.grid(x=0:2, y=0:2)
x <- unlist(mapply(function(x,i) rep(x,i), X$$x, 0:8)) y <- unlist(mapply(function(x,i) rep(x,i), X$$y, 0:8))

Here are the two vectors of $$m=36$$ components:

x 1 2 2 0 0 0 1 1 1 1 2 2 2 2 2 0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2
y 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

This computation provides nice names for the output and prints it:

a <- f(x,y)
rownames(a) <- paste0("x=", 0:2)
colnames(a) <- paste0("y=", 0:2)
a

The output is

y=0 y=1 y=2
x=0   0   3   6
x=1   1   4   7
x=2   2   5   8

As a check, compare this to the output of the table function, table(x,y):

y
x   0 1 2
0 0 3 6
1 1 4 7
2 2 5 8

I was surprised to find that applying table in this instance is an order of magnitude slower than f. I increased $$n$$ until it was equal to $$10001$$ and the difference grew even worse. Apparently, the built-in efficiency of matrix operations in R really makes a difference even for this basic operation!

• This crosstabulation, as I see, amounts to multiplying Gx' Gy where Gx is the dummy set (i.e. design matrix) from variable X and Gy likewise from variable Y. But in your example the values (0,1,2) are known/specified in advance and are same for both variables, therefore, to get Gx and Gy, one don't need to explore the values by tabulation and recoding into dummies (by some function like "design(X)" or whatever it is called). – ttnphns Oct 29 '18 at 16:05
• @ttnphns Yes, but a smart algorithm will scan once over the input vectors at $O(n)$ cost to determine which values to tabulate, so that makes no difference asymptotically. If arbitrary values are to be handled then they will have to be hashed into an associative array at a constant cost of $O(\log n),$ but that cannot wholly explain the consistent inferiority of table compared to f. With a bit of work using table, though, one can achieve good performance: n<-2; matrix(tabulate(1 + x + (n+1)*y, nbins=(n+1)^2), n+1) – whuber Oct 29 '18 at 16:21