Matrix of vector-by-vector sum-squared deviations of two matrices of column vectors

Question

Context

I'm working on a Python program in which I will calculate some number $r$ of matrices $\mathbf{A}^i$ with identical dimensions $m\times n$. For this application, each matrix is probably best thought of as a collection of $L^1$-normed column vectors $\mathbf{c}^i_j$ with strictly non-negative elements,

$$ \mathbf{A}^i = \left[ \begin{array}{ccc} | & \dots & | \\ \mathbf{c}^i_1 & \dots & \mathbf{c}^i_n \\ | & \dots & | \end{array} \right], $$

where each element $c^i_{j,k}$ represents a fractional contribution of one of $m$ components to one of $n$ composite objects. For underlying structural reasons, I anticipate each column vector $\mathbf{c}^i_b$ to have a fairly unambiguous 'closest-match' column vector in both adjacent matrices, which may or may not fall at the same index $b$ in those neighboring matrices:

$$ \mathbf{c}^{i-1}_a \longleftrightarrow \mathbf{c}^i_b \longleftrightarrow \mathbf{c}^{i+1}_c \\ a \overset{?}{=} b \overset{?}{=}c $$

I need to identify the complete set of these 'closest match' mappings.

At present, I am planning to use the element-wise sum-squared deviation (equivalently, since the dimension $m$ is constant among all the matrices, the mean-squared-deviation) between the pairs of column vectors in adjacent matrices as the metric for 'closest match':

$$ \delta^{i}_{b,c} = \sum_k{\left(c^i_{k,b}-c^{i+1}_{k,c}\right)^2} $$

I plan to construct the $r-1$ matrices $\mathbf{\Delta}^i$ of these SSD values,

$$ \mathbf{\Delta}^i = \left[ \begin{array}{ccc} \delta^i_{1,1} & \dots & \delta^i_{1,n} \\ \vdots & \ddots & \vdots \\ \delta^i_{n,1} & \dots & \delta^i_{n,n} \\ \end{array} \right], $$

and then perform further manipulation on the $\mathbf{\Delta}^i$ to identify the desired mappings from each $\mathbf{A}^i$ to the next.

Question

Given such a set of matrices $\mathbf{A}^i$, what is an efficient way to generate the matrices $\mathbf{\Delta}^i$? I would like to exploit basic linear algebra (transpose, matrix multiplication, etc.) to the extent possible. Alternatively, are there other 'closest-match' metrics I might use for which basic linear algebra would be well suited?

Answer 1

The best way I've thought of to achieve this using the sum-squared deviation as the 'closest match' metric starts with rewriting the sum for $\delta^i_{b,c}$ by expanding the polynomial:

$$ \delta^i_{b,c} = \sum_{k}{\left[ \left(c^i_{k,b}\right)^2 - 2c^i_{k,b}c^{i+1}_{k,c} + \left(c^{i+1}_{k,c}\right)^2 \right]} $$

The sums can then be split and reordered,

$$ \delta^i_{b,c} = \sum_{k}{\left(c^i_{k,b}\right)^2} + \sum_{k}{\left(c^{i+1}_{k,c}\right)^2} - 2 \sum_{k}{c^i_{k,b}c^{i+1}_{k,c}}, $$

and the matrix $\mathbf{\Delta}^i$ written as the sum of three contributions, all $n\times n$ matrices:

$$ \mathbf{\Delta}^i = \mathbf{X}^i + \mathbf{Y}^i - 2\mathbf{Z}^i $$

The $\mathbf{Z}^i$ piece is the easiest, since by the definitions of matrix multiplication and transposition:

$$ \mathbf{Z}^i = \left(\mathbf{A}^i\right)^\mathrm{T}\mathbf{A}^{i+1} $$

Dealing with $\mathbf{X}^i$ and $\mathbf{Y}^i$ is problematic, though. The form of both terms is quite familiar: each is the vector dot product of the relevant column vector $\mathbf{c}^i_b$ $\left(\mathbf{c}^{i+1}_c\right)$ with itself. However, each $\mathbf{c}^i_b\cdot\mathbf{c}^i_b \left(\mathbf{c}^{i+1}_c\cdot\mathbf{c}^{i+1}_c\right)$ contributes to every element in row $b$ (column $c$) of $\mathbf{\Delta}^i$.

One way of writing $\mathbf{X}^i$ and $\mathbf{Y}^i$ is as follows:

$$ \mathbf{X}^i = \left[\begin{array}{ccc} \dots & | & \dots \\ \dots & \mathbf{d}^i & \dots \\ \dots & | & \dots \end{array}\right] \\ \mathbf{Y}^i = \left[\begin{array}{ccc} \vdots & \vdots & \vdots \\ - & \mathbf{d}^{i+1} & - \\ \vdots & \vdots & \vdots \end{array}\right] $$

Here, $\mathbf{d}^i$ is the column vector of dot products of the vectors $\mathbf{c}^i_j$, such that $d^i_j=\mathbf{c}^i_j\cdot\mathbf{c}^i_j$. How to obtain these vectors by 'basic' linear algebra, if possible, is not clear to me; nor was it clear to another, as-yet-unanswered questioner two years ago. Fortunately, as I am working in Python, numpy.diagonal will readily extract the vectors $\mathbf{d}^i$ from $\left(\mathbf{A}^i\right)^\mathrm{T}\mathbf{A}^i$. Also fortunately, Python's numpy.repeat is available as a means to stack together replicates of the $\mathbf{d}^i$ vectors into the needed matrices $\mathbf{X}^i$ and, after transposition, $\mathbf{Y}^i$.

While this approach will suit my needs, I'm wondering if there's a better way to obtain $\mathbf{X}^i$ and $\mathbf{Y}^i$, without performing 'matrix surgery' of this sort, to facilitate such calculations in languages with basic linear algebra support but without a rich set of matrix manipulation utilities like Python or MATLAB.