Ferkans — Interactive Telecom Tutor

Why Kronecker Structure Is the Key to Scalable RF Imaging

In Chapter 6 we derived the sensing matrix $\mathbf{A} \in \mathbb{C}^{M \times N}$ that maps the discretized reflectivity $\mathbf{c} \in \mathbb{C}^N$ to measurements $\mathbf{y} \in \mathbb{C}^M$ . For a typical 3D imaging scenario with $N_t = 8$ transmitters, $N_r = 16$ receivers, $N_f = 64$ subcarriers, and a $32^3$ voxel grid, the full matrix has $M \times N = 8192 \times 32768 \approx 2.7 \times 10^8$ complex entries --- over 2 GB of storage. Forming this matrix explicitly, let alone inverting it, is computationally prohibitive.

The central insight of this section is that under physically reasonable assumptions, $\mathbf{A}$ decomposes as a Kronecker product of three much smaller factor matrices. This factorization reduces storage by five orders of magnitude and enables matrix-vector products that are fast enough for iterative reconstruction in real time. Every reconstruction algorithm in Parts IV--VI exploits this structure.

Definition:
Kronecker Product

The Kronecker product of matrices $\mathbf{B} \in \mathbb{C}^{m_1 \times n_1}$ and $\mathbf{C} \in \mathbb{C}^{m_2 \times n_2}$ is the $m_1 m_2 \times n_1 n_2$ block matrix

$\mathbf{B} \otimes \mathbf{C} = \begin{bmatrix} b_{11}\mathbf{C} & b_{12}\mathbf{C} & \cdots & b_{1n_1}\mathbf{C} \\ b_{21}\mathbf{C} & b_{22}\mathbf{C} & \cdots & b_{2n_1}\mathbf{C} \\ \vdots & \vdots & \ddots & \vdots \\ b_{m_1 1}\mathbf{C} & b_{m_1 2}\mathbf{C} & \cdots & b_{m_1 n_1}\mathbf{C} \end{bmatrix}.$

The Kronecker product satisfies the fundamental vec identity:

$(\mathbf{B} \otimes \mathbf{C})\,\text{vec}(\mathbf{X}) = \text{vec}(\mathbf{C}\mathbf{X}\mathbf{B}^T)$

for any conformable matrix $\mathbf{X}$ .

The vec identity is the computational engine behind every fast algorithm in this chapter. It converts a single large matrix-vector product into a sequence of smaller matrix multiplications.

Kronecker product

A matrix operation $\mathbf{B} \otimes \mathbf{C}$ that replaces each element $b_{ij}$ of $\mathbf{B}$ with the block $b_{ij}\mathbf{C}$ , producing an $(m_1 m_2) \times (n_1 n_2)$ matrix from factors of sizes $m_1 \times n_1$ and $m_2 \times n_2$ .

Vectorization operator

The operator $\text{vec}(\mathbf{X})$ that stacks the columns of a matrix $\mathbf{X} \in \mathbb{C}^{m \times n}$ into a single vector in $\mathbb{C}^{mn}$ . The inverse operation reshapes a vector back into a matrix of specified dimensions.

Definition:
Algebraic Properties of the Kronecker Product

The Kronecker product satisfies the following identities, each of which we exploit for computational gain:

Mixed-product rule: $(\mathbf{B}_1 \otimes \mathbf{C}_1)(\mathbf{B}_2 \otimes \mathbf{C}_2) = (\mathbf{B}_1\mathbf{B}_2) \otimes (\mathbf{C}_1\mathbf{C}_2)$

Transpose: $(\mathbf{B} \otimes \mathbf{C})^T = \mathbf{B}^T \otimes \mathbf{C}^T, \qquad (\mathbf{B} \otimes \mathbf{C})^H = \mathbf{B}^H \otimes \mathbf{C}^H$

Inverse (when both factors are invertible): $(\mathbf{B} \otimes \mathbf{C})^{-1} = \mathbf{B}^{-1} \otimes \mathbf{C}^{-1}$

Pseudoinverse: $(\mathbf{B} \otimes \mathbf{C})^\dagger = \mathbf{B}^\dagger \otimes \mathbf{C}^\dagger$

Eigenvalues: If $\mathbf{B}$ has eigenvalues $\{\lambda_i\}$ and $\mathbf{C}$ has eigenvalues $\{\mu_j\}$ , then $\mathbf{B} \otimes \mathbf{C}$ has eigenvalues $\{\lambda_i \mu_j\}$ .

The mixed-product rule is particularly powerful: it means that $\mathbf{A}^{H}\mathbf{A}$ inherits the Kronecker structure of $\mathbf{A}$ , so the Gram matrix can be analyzed factor by factor.

Theorem: Kronecker Factorization of the Sensing Matrix

Consider a multi-static RF imaging system with $N_t$ transmitters, $N_r$ receivers, and $N_f$ frequency bins, operating in the far field with the scene discretized on a separable Cartesian grid of dimensions $N_x \times N_y \times N_r$ . Under these assumptions, the sensing matrix factors as:

$\boxed{\mathbf{A} = \mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}} \otimes \mathbf{A}_{f}}$

where:

$\mathbf{A}_{\text{Rx}} \in \mathbb{C}^{N_r \times N_x}$ : receive spatial factor (columns are Rx steering vectors evaluated at grid angles),
$\mathbf{A}_{\text{Tx}} \in \mathbb{C}^{N_t \times N_y}$ : transmit spatial factor (columns are Tx steering vectors evaluated at grid angles),
$\mathbf{A}_{f} \in \mathbb{C}^{N_f \times N_r}$ : frequency (range) factor (a partial DFT matrix sampling the delay-frequency relationship).

The total matrix $\mathbf{A} \in \mathbb{C}^{M \times N}$ with $M = N_rN_t N_f$ and $N = N_x N_y N_r$ is specified by only $N_r N_x + N_t N_y + N_f N_r$ parameters instead of $MN$ .

The physical reason for separability is that in the far field, the round-trip phase from transmitter $i$ through voxel $(i_x, i_y, i_r)$ to receiver $j$ at frequency $f_k$ separates into three independent terms: one depending only on the Rx-angle pair, one on the Tx-angle pair, and one on the frequency-range pair. Since the sensing matrix entry is the product of these three phase terms, the full matrix is their Kronecker product.

Proof

Separability of the far-field phase

From the Born approximation (Ch 06, Eq. 21), the $(m, n)$ entry of $\mathbf{A}$ for measurement index $m = (j, i, k)$ and voxel $n = (i_x, i_y, i_r)$ involves the phase

$\phi = \kappa_{k}\bigl(d(\mathbf{s}_{i}, \mathbf{p}_{n}) + d(\mathbf{p}_{n}, \mathbf{r}_{j})\bigr).$

For a separable grid $\mathbf{p}_{n} = (x_{i_x}, y_{i_y}, r_{i_r})$ and transmit/receive arrays oriented along orthogonal axes, the far-field Taylor expansion of Ch 06 gives:

$d(\mathbf{s}_{i}, \mathbf{p}_{n}) \approx d(\mathbf{s}_{i}, \mathbf{p}_{0}) + \nabla_{\mathbf{p}} d\big|_{\mathbf{p}_{0}}^T (\mathbf{p}_{n} - \mathbf{p}_{0}).$

The gradient separates into components along the grid axes.

Three independent phase factors

After the Taylor expansion, the phase factors as:

$e^{j\phi} = \underbrace{e^{-j\kappa_{k} \hat{\mathbf{r}}_j^T \mathbf{e}_x \, x_{i_x}}}_{[\mathbf{A}_{\text{Rx}}]_{j, i_x}} \cdot \underbrace{e^{j\kappa_{k} \hat{\mathbf{s}}_i^T \mathbf{e}_y \, y_{i_y}}}_{[\mathbf{A}_{\text{Tx}}]_{i, i_y}} \cdot \underbrace{e^{-j2\pi f_k \tau_{i_r}}}_{[\mathbf{A}_{f}]_{k, i_r}}$

where $\hat{\mathbf{s}}_i$ and $\hat{\mathbf{r}}_j$ are unit direction vectors, and $\tau_{i_r}$ is the round-trip delay for range bin $i_r$ .

Kronecker product identification

Since $[\mathbf{A}]_{(j,i,k),(i_x,i_y,i_r)} = [\mathbf{A}_{\text{Rx}}]_{j,i_x} \cdot [\mathbf{A}_{\text{Tx}}]_{i,i_y} \cdot [\mathbf{A}_{f}]_{k,i_r}$ , this is precisely the definition of the Kronecker product:

$\mathbf{A} = \mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}} \otimes \mathbf{A}_{f}. \quad \blacksquare$

On the Factor Ordering Convention

The ordering $\mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}} \otimes \mathbf{A}_{f}$ corresponds to the measurement index ordering $(j, i, k)$ --- receive antenna varies slowest, frequency varies fastest --- and scene index ordering $(i_x, i_y, i_r)$ . Different stacking conventions lead to permuted Kronecker factors. The CommIT simulator uses the convention above, matching Caire's note (Eq. 22--23). When reading other references, verify the index ordering before comparing expressions.

Definition:
Storage and Computational Savings

For a Kronecker-factored sensing matrix $\mathbf{A} = \mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}} \otimes \mathbf{A}_{f}$ with factor sizes $m_k \times n_k$ :

Quantity	Full $\mathbf{A}$	Kronecker factors
Storage	$M \cdot N$ complex numbers	$\sum_{k=1}^{3} m_k n_k$ complex numbers
Matvec	$O(MN)$ multiplications	$O(m_1 n_1 n_2 n_3 + m_1 m_2 n_2 n_3 + m_1 m_2 m_3 n_3)$

For the representative parameters $N_t = 8$ , $N_r = 16$ , $N_f = 64$ , and $N_x = N_y = N_r = 32$ :

	Full	Kronecker	Ratio
Storage	$2.1$ GB	$18$ KB	$\sim 10^5 \times$
Matvec flops	$2.7 \times 10^8$	$5.2 \times 10^5$	$\sim 500 \times$

The Kronecker factorization makes iterative reconstruction feasible for 3D imaging problems that would otherwise require supercomputer-scale resources.

Theorem: Fast Matrix-Vector Product via Kronecker Structure

For $\mathbf{A} = \mathbf{A}_{3} \otimes \mathbf{A}_{2} \otimes \mathbf{A}_{1}$ with factors of size $m_k \times n_k$ ( $k = 1, 2, 3$ ), the product $\mathbf{y} = \mathbf{A}\mathbf{c}$ can be computed in

$O(m_1 n_1 n_2 n_3 + m_1 m_2 n_2 n_3 + m_1 m_2 m_3 n_3)$

operations, compared to $O(m_1 m_2 m_3 \cdot n_1 n_2 n_3)$ for the naive product. For balanced factors ( $m_k \approx n_k \approx n$ ), this reduces from $O(n^6)$ to $O(n^4)$ .

Instead of applying the giant matrix all at once, we reshape the reflectivity vector into a 3D tensor and apply each factor matrix along one mode at a time --- three small multiplications instead of one enormous one.

Proof

Reshape into a tensor

Reshape $\mathbf{c} \in \mathbb{C}^{n_1 n_2 n_3}$ into a 3D tensor $\boldsymbol{\mathcal{C}} \in \mathbb{C}^{n_1 \times n_2 \times n_3}$ .

Sequential mode products

Apply each factor along one tensor mode:

$\boldsymbol{\mathcal{Y}}_1 = \boldsymbol{\mathcal{C}} \times_1 \mathbf{A}_{1}$ : cost $O(m_1 n_1 n_2 n_3)$ .
$\boldsymbol{\mathcal{Y}}_2 = \boldsymbol{\mathcal{Y}}_1 \times_2 \mathbf{A}_{2}$ : cost $O(m_1 m_2 n_2 n_3)$ .
$\boldsymbol{\mathcal{Y}}_3 = \boldsymbol{\mathcal{Y}}_2 \times_3 \mathbf{A}_{3}$ : cost $O(m_1 m_2 m_3 n_3)$ .

The result $\text{vec}(\boldsymbol{\mathcal{Y}}_3) = \mathbf{A}\mathbf{c}$ . $\blacksquare$

Example: Kronecker Structure for 2D Imaging with ULAs

A monostatic system uses a $N_t = 4$ -element Tx ULA and a co-located $N_r = 8$ -element Rx ULA, both with half-wavelength spacing, transmitting on $N_f = 16$ subcarriers with spacing $\Delta f = 1$ MHz centered at $f_0 = 28$ GHz. The scene is a 2D grid of $N_x = 16$ cross-range bins and $N_r = 16$ range bins.

(a) Write the factor matrices explicitly. (b) Compute the storage saving. (c) Determine the computational saving for one matvec.

Solution

Factor matrices

Since this is a 2D problem, $\mathbf{A} = \mathbf{A}_{\text{ang}} \otimes \mathbf{A}_{f}$ where the angular factor combines Tx and Rx via the virtual array (Ch 06.7). The virtual array has $N_tN_r = 32$ elements.

$[\mathbf{A}_{\text{ang}}]_{v, i_x} = e^{-j\kappa d_v \sin\theta_{i_x}}, \quad v = 1,\ldots,32, \; i_x = 1,\ldots,16$

$[\mathbf{A}_{f}]_{k, i_r} = e^{-j2\pi f_k \tau_{i_r}}, \quad k = 1,\ldots,16, \; i_r = 1,\ldots,16$

So $\mathbf{A}_{\text{ang}} \in \mathbb{C}^{32 \times 16}$ and $\mathbf{A}_{f} \in \mathbb{C}^{16 \times 16}$ .

Storage

Full matrix: $M \times N = 512 \times 256 = 131{,}072$ complex entries $= 2$ MB (at 16 bytes per complex double).

Kronecker factors: $32 \times 16 + 16 \times 16 = 768$ entries $= 12$ KB. Saving: $\mathbf{170\times}$ .

Matvec cost

Full matvec: $2 \times 512 \times 256 = 262{,}144$ flops.

Kronecker matvec: reshape $\mathbf{c}$ into $16 \times 16$ , then $\mathbf{A}_{f} \mathbf{X}$ costs $2 \times 16 \times 16 \times 16 = 8192$ flops, $(\mathbf{A}_{f} \mathbf{X})\mathbf{A}_{\text{ang}}^T$ costs $2 \times 16 \times 32 \times 16 = 16{,}384$ flops. Total: $24{,}576$ flops. Saving: $\mathbf{10.7\times}$ .

When the factors are DFT matrices (uniform grid, uniform frequency spacing), each product becomes an FFT: $O(n \log n)$ per mode.

Kronecker Structure Visualization

Visualizes the magnitude of the sensing matrix and its Kronecker factors.

Top row: Magnitude plots of the three factor matrices $\mathbf{A}_{\text{Rx}}$ , $\mathbf{A}_{\text{Tx}}$ , $\mathbf{A}_{f}$ .

Bottom left: Magnitude of the full sensing matrix $\mathbf{A}$ , showing the characteristic block-repetitive pattern.

Bottom right: Singular values of $\mathbf{A}$ (blue) compared with products of factor singular values (red circles), confirming the multiplicative SVD property.

Parameters

N_t

(Tx antennas)4

N_r

(Rx antennas)8

N_f

(Frequencies)16

N_x

(Grid size per dimension)12

The Gram Matrix Inherits Kronecker Structure

Since $\mathbf{A} = \mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}} \otimes \mathbf{A}_{f}$ , the Gram matrix factors as

$\mathbf{A}^{H}\mathbf{A} = (\mathbf{A}_{\text{Rx}}^H\mathbf{A}_{\text{Rx}}) \otimes (\mathbf{A}_{\text{Tx}}^H\mathbf{A}_{\text{Tx}}) \otimes (\mathbf{A}_{f}^{H}\mathbf{A}_{f}).$

This means the point-spread function (PSF) --- the response of the imaging system to a point scatterer --- can be analyzed factor by factor. The overall PSF is the Kronecker product of the PSFs in the receive-angle, transmit-angle, and range dimensions. A narrow main lobe in one dimension cannot compensate for wide sidelobes in another: the imaging resolution is determined by the worst factor.

Common Mistake: When Kronecker Structure Breaks Down

Mistake:

Assuming exact Kronecker decomposition in all operating conditions. The factorization requires:

Far-field approximation --- spherical wavefronts in the near field break separability.
Separable Cartesian grid --- polar or adaptive grids prevent clean factorization.
Narrowband per subcarrier --- wideband effects within a single subcarrier couple range and angle.
No mutual coupling --- antenna interactions create non-separable effects.

Correction:

When these assumptions are only approximately satisfied, $\mathbf{A} \approx \mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}} \otimes \mathbf{A}_{f} + \mathbf{E}$ where $\mathbf{E}$ is a perturbation. Low-rank approximations of $\mathbf{E}$ or iterative refinement can recover most of the computational benefit. For near-field scenarios, see the extended Fresnel-zone corrections in Ch 08.4.

Historical Note: Leopold Kronecker and the Product That Bears His Name

1858--present

The Kronecker product was introduced by the German mathematician Leopold Kronecker (1823--1891), though it was also independently developed by Johann Georg Zehfuss in 1858. The operation was initially a curiosity in pure algebra. Its computational significance was recognized much later, when Van Loan and Pitsianis (1993) showed that Kronecker product approximation provides optimal low-rank factorizations for structured matrices --- precisely the setting we encounter in imaging. In the radar and signal processing community, the Kronecker structure was exploited for MIMO radar by Li and Stoica (2007), and Caire's framework extends it to the multi-static, multi-frequency RF imaging context.

Separable grid

A discretization of the target region where the voxel positions form a Cartesian product: $\{(x_i, y_j, r_k) : i=1,\ldots,N_x;\; j=1,\ldots,N_y;\; k=1,\ldots,N_r\}$ . This separability is what enables the Kronecker factorization of $\mathbf{A}$ .

🔧Engineering Note

GPU Implementation of Kronecker Matvec

The sequential mode-product algorithm maps naturally to GPU computation. Each mode product is a batched matrix multiplication:

Mode 1 ( $\mathbf{A}_{f}$ along range): batch of $N_x N_y$ matrix-vector products of size $N_f \times N_r$ .
Mode 2 ( $\mathbf{A}_{\text{Tx}}$ along Tx-angle): batch of $N_x m_1$ products of size $N_t \times N_y$ .
Mode 3 ( $\mathbf{A}_{\text{Rx}}$ along Rx-angle): batch of $m_1 m_2$ products of size $N_r \times N_x$ .

Using CuPy or PyTorch, each mode product is a single torch.einsum or cupy.tensordot call, achieving near-peak GPU throughput. The CommIT simulator implements this pattern, enabling real-time imaging at millimeter-wave frequencies.

Key Takeaway

Under far-field and separable-grid assumptions, the sensing matrix factors as $\mathbf{A} = \mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}} \otimes \mathbf{A}_{f}$ . This reduces storage from $O(MN)$ to $O(\sum m_k n_k)$ (a factor of $\sim 10^5$ for typical parameters) and matrix-vector products from $O(MN)$ to sequential mode products costing $O(N^{4/3})$ for balanced factors, or $O(N \log N)$ when factors are DFT matrices. Every reconstruction algorithm in this book exploits this structure.

Kronecker Factorization of the Sensing Matrix

Animated decomposition of the full sensing matrix

\mathbf{A}

into its Kronecker factors

\mathbf{A}_{\text{freq}} \otimes \mathbf{A}_{\text{Rx}} \otimes \mathbf{A}_{\text{Tx}}

. The animation shows how the factor dimensions relate to the physical array and frequency grid, and how matrix-vector products factorize into a sequence of smaller operations — the key insight behind efficient RF imaging algorithms.

Kronecker Product Structure of the Sensing Matrix