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? 
 A: 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)$$

is the answer.
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!
