---
title: Khatri-Rao Random Projections (KRPs)
url: https://www.emergentmind.com/topics/khatri-rao-random-projections-krps
type: topic
---

# Khatri-Rao Random Projections (KRPs)

Searching arXiv for recent and foundational papers on Khatri-Rao random projections.
Khatri–Rao Random Projections (KRPs) are structured random embeddings built from Khatri–Rao products of smaller random matrices. In this construction, the sketching matrix is assembled columnwise as a Kronecker product of corresponding columns from multiple factors, yielding a map that preserves key geometric or inverse-problem structure while reducing storage and sketch-generation cost relative to dense unstructured random matrices. Across compressed sensing, Johnson–Lindenstrauss-type dimension reduction, randomized low-rank approximation, tensor compression, eigensolvers, and tensor-train rounding, the central theme is that Khatri–Rao structure can provide either stronger restricted isometry than the constituent factors or substantially lower memory and arithmetic cost when the ambient problem already exhibits Kronecker or tensor-product structure [1709.05789], [2105.00105], [2507.23207], [2405.11962].

## 1. Definition and basic construction

For matrices \(A=[a_1,\dots,a_n]\in\mathbb{R}^{m\times n}\) and \(B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}\), their columnwise Khatri–Rao product \(A\odot B\in\mathbb{R}^{m^2\times n}\) is defined by
\[
(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,
\]
where \(\otimes\) denotes the Kronecker product [1709.05789]. More generally, if \(A\in\mathbb{R}^{I\times K}\) and \(B\in\mathbb{R}^{J\times K}\), then \(A\odot B\in\mathbb{R}^{(IJ)\times K}\) has \(\ell\)-th column \(A(:,\ell)\otimes B(:,\ell)\), and the operation is associative, so \((A\odot B)\odot C=A\odot(B\odot C)\) [2105.00105]. For \(d\) factor matrices \(\Omega^{(j)}\in\mathbb{R}^{n_j\times \ell}\), \(j=1,\dots,d\), one obtains the Khatri–Rao sketch
\[
\Omega=\Omega^{(1)}\odot\Omega^{(2)}\odot\cdots\odot\Omega^{(d)}\in\mathbb{R}^{N\times \ell},\qquad N=\prod_j n_j,
\]
whose columns are isotropic random vectors when the factors have i.i.d. mean-zero, unit-variance subgaussian entries [2507.23207].

In dimension-reduction form, the embedding is linear:
\[
f_{\rm TRP}(x)=\Phi^T x\in\mathbb{R}^k,
\]
where \(\Phi=A_1\odot\cdots\odot A_N\) is formed from smaller random projections \(A_n\in\mathbb{R}^{d_n\times k}\) satisfying \(d=\prod_{n=1}^N d_n\) [2105.00105]. In the tensor viewpoint, each row of \(\Phi\) corresponds to a multi-index in the tensorized ambient space. This factorized construction is the basis of Tensor Random Projection (TRP), which is one of the main KRP-based random projection models for low-memory embedding [2105.00105].

A recurrent algebraic identity is the Rao–Rao relation
\[
(A_S\odot B_S)^T(A_S\odot B_S)=(A_S^T A_S)\circ(B_S^T B_S),
\]
with \(\circ\) the Hadamard product [1709.05789]. This identity explains why the geometry of KRP sketches is governed by Hadamard products of correlation or Gram matrices, and it underlies both restricted isometry analyses and algorithmic accelerations exploiting Kronecker structure.

## 2. Restricted isometry and sparse recovery

One of the earliest systematic analyses of KRPs in compressed sensing studies the restricted isometry property (RIP) of columnwise Khatri–Rao products [1709.05789]. For a matrix \(\Phi\in\mathbb{R}^{M\times N}\), the \(k\)-restricted isometry constant \(\delta_k(\Phi)\in[0,1)\) is defined by
\[
(1-\delta_k)\|z\|_2^2\le \|\Phi z\|_2^2\le (1+\delta_k)\|z\|_2^2
\]
for every \(k\)-sparse \(z\in\mathbb{R}^N\), with the smallest such \(\delta_k\) called the \(k\)-RIC [1709.05789].

A deterministic result establishes that if \(A,B\in\mathbb{R}^{m\times n}\) have unit-\(\ell_2\)-norm columns and \(k\)-RICs \(\delta_k(A),\delta_k(B)<1\), and if \(\delta=\max\{\delta_k(A),\delta_k(B)\}\), then
\[
\delta_k(A\odot B)\le \delta^2.
\]
Equivalently, for every \(k\)-sparse \(z\),
\[
(1-\delta^2)\|z\|_2^2\le \|(A\odot B)z\|_2^2\le (1+\delta^2)\|z\|_2^2.
\]
Thus the Khatri–Rao product exhibits stronger restricted isometry than either constituent matrix at the same sparsity level [1709.05789]. The proof proceeds by bounding eigenvalues of the Hadamard product Gram matrix through forward and reverse matrix Kantorovich inequalities.

In the subgaussian setting, if \(A,B\in\mathbb{R}^{m\times n}\) have i.i.d. entries with mean zero, variance one, and subgaussian \(\psi_2\)-norm at most \(\kappa\), and
\[
\Phi=(A/\sqrt{m})\odot(B/\sqrt{m}),
\]
then there exist universal constants \(c,C>0\) depending on \(\kappa\) such that, for fixed \(\delta\in(0,1)\) and \(\gamma>1\), the condition
\[
m\ge C\cdot \gamma\cdot \kappa^4\cdot \frac{k\log n}{\delta}
\]
implies \(\delta_k(\Phi)\le \delta\) with probability at least \(1-10n^{-2(\gamma-1)}\) [1709.05789]. An analogous self-Khatri–Rao result for \(\Phi=(A/\sqrt{m})\odot(A/\sqrt{m})\) holds with probability at least \(1-5n^{-2(\gamma-1)}\) under the same scaling [1709.05789]. Quantitatively, the constituent matrix \(A/\sqrt{m}\) satisfies \(\delta_k(A/\sqrt{m})=O(\sqrt{(k\log n)/m})\), whereas the Khatri–Rao product satisfies \(\delta_k(\Phi)=O((k\log n)/m)\). This yields \(m=O((k\log n)/\delta)\) rather than \(m=O((k\log n)/\delta^2)\) for attaining \(\delta_k\le \delta\) [1709.05789].

The centered self Khatri–Rao setting sharpens the picture for quadratic or covariance-type sensing models. If \(A=[a_1,\dots,a_N]\in\mathbb{R}^{n\times N}\) has iid columns drawn uniformly from the sphere of radius \(\sqrt{n}\) or with iid zero-mean, unit-variance sub-Gaussian entries, the self product has columns \(a_i\otimes a_i=\mathrm{vec}(a_i a_i^T)\). Because \(E[a_i\otimes a_i]=\mathrm{vec}(I_n)\), one defines the centered matrix \(\mathcal{A}\) with columns
\[
\mathcal{A}_i=\kappa(n)\bigl(a_i\otimes a_i-\mathrm{vec}(I_n)\bigr),
\]
and the normalized operator \(\Phi=(1/n)\mathcal{A}\) [1905.09245]. For this centered self Khatri–Rao product, there exist universal constants \(C,c>0\) such that if
\[
s\lesssim \frac{n^2}{\log^2(eN/n^2)},
\]
then
\[
P\{\delta_s(\Phi)\ge \delta\}\le C\exp(-cn)
\]
for any fixed \(\delta\in(0,1)\) [1905.09245]. The analysis uses a heavy-tailed-column RIP theorem together with Hanson–Wright control of the centered quadratic forms. The resulting scaling is applicable to covariance matching, activity detection, and MIMO gain estimation, where measurements are linear in outer products \(a_i a_i^T\) rather than in the vectors \(a_i\) themselves [1905.09245].

These RIP results substantiate a common interpretation of KRPs as “isometry-amplifying” structured sensing matrices. This suggests that, in sparse inverse problems naturally producing second-order or tensorized measurements, KRP structure is not merely a storage device but a geometric advantage.

## 3. Tensor Random Projection and low-memory embeddings

In randomized dimension reduction, Khatri–Rao structure appears as Tensor Random Projection, which constructs a sketch from several smaller random maps instead of a full \(k\times d\) dense matrix [2105.00105]. If \(d=\prod_{n=1}^N d_n\), choose \(A_n\in\mathbb{R}^{d_n\times k}\) and set
\[
\Phi=A_1\odot A_2\odot\cdots\odot A_N\in\mathbb{R}^{(d_1\cdots d_N)\times k}.
\]
The memory cost is then \(O(k\sum_n d_n)\), compared with \(O(kd)\) for a classical dense random projection and \(O(kd\delta)\) for a sparse projection with density \(\delta\) [2105.00105]. If \(d_n\approx d^{1/N}\), this becomes \(O(Nkd^{1/N})\), and for \(N=2\) with \(d=d_1d_2\), storage is approximately \(2k\sqrt{d}\) [2105.00105].

The basic isometry statement is unbiasedness: if each factor has independent columns in isotropic position, then
\[
\mathbb{E}\|f_{\rm TRP}(x)\|^2=\|x\|^2
\]
for every \(x\) [2105.00105]. When the entries of each factor are iid with mean \(0\), variance \(1\), and fourth moment \(\Delta\), the variance of the squared embedding norm is
\[
\mathrm{Var}\bigl(\|f_{\rm TRP}(x)\|^2\bigr)=\frac{1}{k}\Bigl[(\Delta^N-3)\|x\|_4^4+2\|x\|_2^4\Bigr].
\]
A variance-reduced extension, denoted TRP(T), averages \(T\) independent TRP maps and preserves the mean while reducing the first variance term to
\[
\frac{1}{Tk}(\Delta^N-3)\|x\|_4^4+\frac{2}{k}\|x\|_2^4
\]
[2105.00105].

For \(N=2\), a non-asymptotic error bound is available under sub-Gaussian assumptions on the data vector \(x\) and the factor entries:
\[
\mathbb{P}\Bigl(\bigl|\|f_{\rm TRP}(x)\|^2-\|x\|^2\bigr|\ge \varepsilon \|x\|^2\Bigr)\le C\exp\!\bigl[-c_1(\sqrt{k}\,\varepsilon)^{1/2}\bigr].
\]
Via a standard JL-style argument, this yields an \(\varepsilon\)-Johnson–Lindenstrauss embedding for \(n\) points when
\[
k=\mathcal{O}(\varepsilon^{-2}\log^4 n),
\]
whereas dense Gaussian or fast JL methods require \(k=O(\varepsilon^{-2}\log n)\) [2105.00105]. The paper explicitly notes that the \(\log^4 n\) factor is pessimistic empirically.

The empirical results reported for synthetic Gaussian vectors and MNIST emphasize the storage–accuracy trade-off rather than asymptotic optimality. TRP achieves almost the same distortion as its base random projection while using only \(O(2k\sqrt{d})\) memory instead of \(kd\), and TRP(5) nearly matches the variance of dense random projections at approximately \(1/20\) the storage for \(d=40000\) [2105.00105]. Sparse and very sparse choices of the factors preserve the “database-friendly” property of few nonzero queries and no floating-point operations [2105.00105].

A common misconception is that KRPs are only relevant for explicitly tensorized data. The TRP formulation shows a broader principle: tensorization may be imposed on the sketch rather than on the data representation itself, provided the ambient dimension factors as \(\prod_n d_n\). The penalty is looser worst-case embedding bounds; the benefit is a large reduction in memory and random-bit generation.

## 4. Subspace embeddings and low-rank approximation

Recent work has developed KRP guarantees beyond pointwise norm preservation and sparse recovery, notably for oblivious subspace embeddings (OSEs) and randomized range finding [2507.23207], [2405.11962]. These results are central in large-scale matrix and tensor approximation, where the sketch is used to identify a near-optimal subspace rather than to embed individual vectors alone.

For a fixed \(X\in\mathbb{R}^{d\times n}\), a random \(\Omega\in\mathbb{R}^{n\times \ell}\) is an \((\varepsilon,\delta,d)\)-OSE if, with probability at least \(1-\delta\),
\[
(1-\varepsilon)\|Xv\|_2^2\le \|\Omega^T Xv\|_2^2\le (1+\varepsilon)\|Xv\|_2^2
\quad \forall v\in\mathbb{R}^n
\]
[2405.11962]. For the two-factor Gaussian KRP
\[
\Omega=\frac{1}{\sqrt{\ell}}(\tilde\Omega\odot \hat\Omega)\in\mathbb{R}^{n\times \ell},
\qquad n=n_1n_2,
\]
with iid \(N(0,1)\) entries in the factors, one theorem states that \(\Omega\) is an \((\varepsilon,\delta,k)\)-OSE provided
\[
\ell \ge C\Bigl(k^{3/2}\varepsilon^{-2}+k\log(1/\delta)\varepsilon^{-2}
+k^{1/2}\log^2(1/\delta)\varepsilon^{-1}\Bigr),
\qquad C=(2000e^4)^2
\]
[2405.11962]. The dependence on \(k\) is therefore worse than the linear dependence of unstructured Gaussian sketches, but still polynomial rather than catastrophic.

A more general \(d\)-factor subgaussian theory refines this picture. Let \(W\in\mathbb{R}^{N\times r}\) have orthonormal columns and \(\Omega=\Omega^{(1)}\odot\cdots\odot\Omega^{(d)}\in\mathbb{R}^{N\times \ell}\) be a KRP with i.i.d. mean-zero, unit-variance subgaussian entries in each factor. If
\[
\ell \ge 2.6\,\varepsilon^{-2}\,[\,r+rC_{K,d}\ln^d(8\ell/\delta)\,]\ln(4r/\delta),
\]
then with probability at least \(1-\delta\),
\[
(1-\varepsilon)\ell \le \sigma_r^2(\Omega^TW)\le \sigma_1^2(\Omega^TW)\le (1+\varepsilon)\ell,
\]
so \(\Omega/\sqrt{\ell}\) is a \((2\varepsilon)\)-subspace embedding for \(\mathrm{range}(W)\) [2507.23207]. The same paper proves a Hanson–Wright-type tail bound for a single KRP column and a randomized range-finder guarantee: if \(M\in\mathbb{R}^{M\times N}\) has rank-\(r\) leading singular space and the sketch dimension satisfies
\[
\ell \ge 8\,(r+rC_{K,d}\ln^d(8\ell/\delta))\ln(4r/\delta),
\]
then \(Q=\mathrm{orth}(M\Omega)\) obeys
\[
\|M-QQ^TM\|_F^2 \le (1+2r(1+2\Gamma))\|\Sigma_\perp\|_F^2
\]
with explicit \(\Gamma\) depending on \(C_{K,d}\), \(\ell\), \(\delta\), and \(N-r\) [2507.23207].

The significance of these bounds lies less in their raw constants than in the structural improvement over earlier pessimistic analyses. The 2025 analysis emphasizes removal of large explicit factors in \(N\) and only poly-logarithmic dependence on the tensor order \(d\), rather than doubly exponential dependence [2507.23207]. This suggests that much of the previous theoretical gap between empirical and worst-case KRP performance was due to analysis rather than intrinsic instability.

## 5. Structured algorithms in matrix, tensor, and eigenvalue computations

KRP sketches are especially useful when the target matrix or operator already has block, Kronecker, or tensor-product structure. In those regimes, the cost of applying a KRP can be reduced to structured multiplies, Hadamard products, MTTKRP kernels, or partial contractions rather than dense matrix products.

For block-structured matrices \(M\in\mathbb{R}^{mp\times nq}\) of the form
\[
M=\sum_{j=1}^t E_j\otimes M_j,
\]
a single-view randomized SVD using KRPs draws
\[
\Omega=\Omega^{(1)}\odot\Omega^{(2)}\in\mathbb{R}^{nq\times \ell_r},
\qquad
\Psi=\Psi^{(1)}\odot\Psi^{(2)}\in\mathbb{R}^{mp\times \ell_l},
\]
computes structured sketches
\[
Y=M\Omega=\sum_j (E_j\Omega^{(1)})(M_j\Omega^{(2)}),\qquad
Z=\Psi^TM=\sum_j (\Psi^{(1)T}E_j)(\Psi^{(2)T}M_j),
\]
orthogonalizes \(Y=QR\), and recovers a core matrix via \(W=(\Psi^TQ)^\dagger Z\) [2507.23207]. Assuming \(\ell_r,\ell_l=O(r)\), the sketching cost is
\[
O\Bigl(\sum_j r(\mathrm{nnz}(E_j)+\mathrm{nnz}(M_j))\Bigr),
\]
with orthogonalization and least squares costing \(O((mp+nq)r^2)\), compared with \(O(nq\cdot mp\cdot r)\) for dense Gaussian sketches [2507.23207].

For Tucker compression, KRPs replace Gaussian sketches in RHOSVD and sequentially truncated HOSVD. In RHOSVD-KRP, each mode-\(i\) unfolding is multiplied by a KRP over all non-\(i\) modes,
\[
Y_{(i)}=X_{(i)}(\Omega^{(d)}\odot\cdots\odot\Omega^{(i+1)}\odot\Omega^{(i-1)}\odot\cdots\odot\Omega^{(1)}),
\]
followed by thin QR to obtain \(Q_i\), after which the Tucker core is formed by multilinear projection [2507.23207]. This naturally translates sketching into MTTKRP operations. Reported complexities are
\[
O(d\,n\,\ell^2+d\,n^d\,\ell+\sum_{i=1}^d \ell^i n^{d-i+1})
\]
for RHOSVD-KRP and
\[
O(d\,n\,\ell^2+\sum_{i=1}^d 2\ell^i n^{d-i+1})
\]
for RSTHOSVD-KRP under \(n_i\approx n\) and \(\ell\approx r+p\), with memoized RHOSVD-KRP further reducing MTTKRP cost by approximately \(d/2\) [2507.23207].

In eigensolvers for Kronecker-structured operators, KRP sketches serve as structured starting blocks or projection matrices. For contour-integral methods, using \(\Omega=\tilde\Omega\odot\hat\Omega\) enables each projected resolvent application to be reshaped into a Sylvester equation with low-rank right-hand side, which can then be handled by low-rank BiCGstab together with ADI [2405.11962]. In low-rank LOBPCG, if \(A=\sum_{m=1}^s \tilde A_m\otimes \hat A_m\), matrix applications and preconditioner applications can be carried out in block-low-rank tensor format at cost \(O(s(n_1^2r+n_2^2r))\) rather than \(O(n_1^2n_2^2)\) [2405.11962]. The paper explicitly notes that LOBPCG is easier to integrate with KRP structure algorithmically, though with less theoretical justification than the contour-integral setting [2405.11962].

These methods reflect a key distinction between KRP theory and KRP practice. The sketch dimension may be mildly worse than Gaussian in worst-case OSE theory, but the application cost can be drastically lower when the operator structure aligns with the Khatri–Rao factorization.

## 6. Tensor-train rounding, adaptive sketching, and applications

A further development places KRPs at the center of randomized tensor-train (TT) rounding [2511.03598]. In TT format, a \(d\)-way tensor \(\mathcal{X}\) is represented by cores
\[
\mathcal{X}_k\in\mathbb{R}^{r_{k-1}\times n_k\times r_k},\qquad r_0=r_d=1.
\]
Rounding compresses these ranks, traditionally by deterministic TT-SVD or related procedures. The KRP-based approach instead sketches the large unfolding matrices using products
\[
\Omega_{k+1}\odot\Omega_{k+2}\odot\cdots\odot\Omega_d
\in \mathbb{R}^{(\prod_{i=k+1}^d n_i)\times \ell},
\]
where each \(\Omega_i\in\mathbb{R}^{n_i\times \ell}\) is Gaussian [2511.03598].

The practical mechanism is a sequence of right-to-left partial contractions. For \(k=d\),
\[
W_d=H(\mathcal{X}_d)\Omega_d\in\mathbb{R}^{r_{d-1}\times \ell},
\]
and for \(k=d-1,\dots,2\),
\[
W_k=H(\mathcal{X}_k)(W_{k+1}\odot \Omega_k)\in\mathbb{R}^{r_{k-1}\times \ell},
\]
where \(H(\mathcal{X}_k)\) is the horizontal unfolding of the \(k\)-th core [2511.03598]. Each step is an MTTKRP costing \(O(nr^2\ell)\), giving total cost
\[
O(dnr^2\ell)
\]
for the partial contractions and the same leading term overall for fixed-rank KRP-based TT-rounding [2511.03598].

The adaptive variant seeks \(\|\mathcal{X}-\mathcal{Y}\|\le \varepsilon\|\mathcal{X}\|\) without prespecified TT-ranks. The procedure increases sketch block size incrementally, estimates residual norms from KRP sketches, and stops once the local residual drops below
\[
\tau=\varepsilon\|\mathcal{X}\|/\sqrt{d-1}
\]
[2511.03598]. The theoretical basis is a Frobenius-norm estimator: for
\[
\Omega=\Omega_d\odot\cdots\odot\Omega_1\in\mathbb{R}^{N\times \ell},
\]
the quantity
\[
\widehat{\|A\|}_F^2=\frac{1}{\ell}\|A\Omega\|_F^2
\]
is an unbiased estimator of \(\|A\|_F^2\), and if
\[
\ell\ge 4\max\Bigl(\varepsilon^{-1}C''\ln^d\frac{2}{\delta},\;
\varepsilon^{-2}(C')^2\ln^2\frac{2}{\delta}\Bigr),
\]
then with probability at least \(1-\delta\),
\[
\Bigl|\frac{1}{\ell}\|A\Omega\|_F^2-\|A\|_F^2\Bigr|
\le \varepsilon \|A\|_F^2
\]
[2511.03598]. By summing the local errors in quadrature, the final TT approximation satisfies
\[
\|\mathcal{X}-\mathcal{Y}\|_F\le \varepsilon\|\mathcal{X}\|_F
\]
[2511.03598].

The applications reported for KRP-based TT rounding include synthetic low-rank tensors, parametric Matérn kernel approximation, and TT-GMRES for a parametric PDE. The numerical results show speed-ups of up to \(50\times\) over deterministic TT-rounding while maintaining comparable accuracy and only low overhead from adaptivity [2511.03598]. Since the manuscript postdates some earlier KRP theory, this suggests a broadening of the KRP paradigm from static sketching to adaptive stopping and online error control in hierarchical tensor formats.

## 7. Empirical performance, limitations, and research directions

Across the cited literature, numerical evidence consistently shows that KRP sketches are substantially cheaper to generate and apply than dense Gaussian sketches when structure is exploitable, while empirical approximation quality is often similar. In block-Hankel system identification, KRP-based randomized SVD uses approximately \(1.3\%\) of the random numbers of dense Gaussian sketching, runs approximately \(4.5\times\) faster than dense Gaussian sketching and approximately \(3\times\) faster than a block-Hankel-aware Gaussian method, at equal or better accuracy [2507.23207]. In synthetic four-way Cauchy tensor compression, RHOSVD-KRP is approximately \(3\)–\(5\times\) faster than RHOSVD at comparable approximation error, and the memoized version is the fastest among the tested methods [2507.23207]. In tensor-based sensor placement for fluid flow, RHOSVD-KRP-MEMO is approximately \(16\times\) faster than HOSVD and RHOSVD-KRP approximately \(14\times\) faster, at essentially identical test error to RHOSVD [2507.23207]. In eigensolvers, low-rank BiCGstab within a contour-integral framework yields large memory savings relative to direct sparse solves, and low-rank LOBPCG reaches eigenvalue errors below \(10^{-8}\) on \(3000^2\) grids with per-iteration times reported between approximately \(0.45\) s and \(2.0\) s depending on the potential [2405.11962].

The literature also states several limitations. In TRP, variance scales with \(\Delta^N\), so using too many factors degrades accuracy [2105.00105]. The worst-case embedding dimension for TRP has an extra \(\log^3 n\) factor in the \(N=2\) analysis, and the KRP OSE bound in eigensolver-oriented theory depends on \(k^{3/2}\) rather than linearly on \(k\) [2105.00105], [2405.11962]. In eigenvalue computations, algorithmic gains in contour-integral methods depend strongly on the ability to solve shifted matrix equations efficiently and accurately in low-rank form [2405.11962]. For TT rounding, the theoretical guarantees assume Gaussian factors and independence across cores [2511.03598].

Several open directions are explicitly identified. These include deterministic RIP bounds when the factor matrices are not unit-norm or do not themselves satisfy RIP; extension to Khatri–Rao products of more than two matrices with tighter concentration for heavy-tailed entries; and exploitation of Khatri–Rao structure in fast algorithms for large-scale sparse recovery and dimensionality reduction [1709.05789]. Later work effectively pursues these directions by improving dependence on tensor order, integrating KRPs into randomized Tucker algorithms, and developing adaptive error-estimation schemes for TT rounding [2507.23207], [2511.03598].

A plausible synthesis is that KRP research has evolved along two coupled lines. The first concerns geometry: stronger RIP, unbiased norm preservation, and OSE guarantees for structured random embeddings. The second concerns computation: replacing dense sketch formation and application by MTTKRP, Hadamard-Gram identities, Kronecker-aware multiplies, partial contractions, and low-rank Sylvester solves. The current state of the field indicates that KRPs are most effective not as generic replacements for Gaussian sketches in all settings, but as structure-aligned random projections whose advantages become pronounced when the ambient problem already contains tensor, covariance, or Kronecker-product organization [1709.05789], [2105.00105], [2405.11962], [2507.23207], [2511.03598].

Source: https://www.emergentmind.com/topics/khatri-rao-random-projections-krps