Khatri-Rao Random Projections (KRPs)

• Khatri-Rao Random Projections (KRPs) are structured random embeddings constructed via the Khatri–Rao product of smaller matrices, enabling efficient dimension reduction.
• They enhance restricted isometry and promote lower storage costs by leveraging inherent Kronecker or tensor-product structures, making them ideal for compressed sensing and low-rank approximations.
• KRPs underpin fast algorithms in matrix, tensor, and eigenvalue computations by exploiting structure to reduce arithmetic operations and memory footprint.

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 (Khanna et al., 2017, Sun et al., 2021, Saibaba et al., 31 Jul 2025, Bujanović et al., 2024).

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 (Khanna et al., 2017). 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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$0 (Sun et al., 2021). For $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$1 factor matrices $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$2, $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$3, one obtains the Khatri–Rao sketch

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$4

whose columns are isotropic random vectors when the factors have i.i.d. mean-zero, unit-variance subgaussian entries (Saibaba et al., 31 Jul 2025).

In dimension-reduction form, the embedding is linear: $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$5 where $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$6 is formed from smaller random projections $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$7 satisfying $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$8 (Sun et al., 2021). In the tensor viewpoint, each row of $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$9 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 (Sun et al., 2021).

A recurrent algebraic identity is the Rao–Rao relation

$A\odot B\in\mathbb{R}^{m^2\times n}$0

with $A\odot B\in\mathbb{R}^{m^2\times n}$1 the Hadamard product (Khanna et al., 2017). 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 (Khanna et al., 2017). For a matrix $A\odot B\in\mathbb{R}^{m^2\times n}$2, the $A\odot B\in\mathbb{R}^{m^2\times n}$3-restricted isometry constant $A\odot B\in\mathbb{R}^{m^2\times n}$4 is defined by

$A\odot B\in\mathbb{R}^{m^2\times n}$5

for every $A\odot B\in\mathbb{R}^{m^2\times n}$6-sparse $A\odot B\in\mathbb{R}^{m^2\times n}$7, with the smallest such $A\odot B\in\mathbb{R}^{m^2\times n}$8 called the $A\odot B\in\mathbb{R}^{m^2\times n}$9-RIC (Khanna et al., 2017).

A deterministic result establishes that if $$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$0 have unit-$$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$1-norm columns and $$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$2-RICs $$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$3, and if $$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$4, then

$$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$5

Equivalently, for every $$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$6-sparse $$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$7,

$$(A\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$8

Thus the Khatri–Rao product exhibits stronger restricted isometry than either constituent matrix at the same sparsity level (Khanna et al., 2017). 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\odot B)_{:,i}=a_i\otimes b_i,\qquad i=1,\dots,n,$$9 have i.i.d. entries with mean zero, variance one, and subgaussian $\otimes$0-norm at most $\otimes$1, and

$\otimes$2

then there exist universal constants $\otimes$3 depending on $\otimes$4 such that, for fixed $\otimes$5 and $\otimes$6, the condition

$\otimes$7

implies $\otimes$8 with probability at least $\otimes$9 (Khanna et al., 2017). An analogous self-Khatri–Rao result for $A\in\mathbb{R}^{I\times K}$0 holds with probability at least $A\in\mathbb{R}^{I\times K}$1 under the same scaling (Khanna et al., 2017). Quantitatively, the constituent matrix $A\in\mathbb{R}^{I\times K}$2 satisfies $A\in\mathbb{R}^{I\times K}$3, whereas the Khatri–Rao product satisfies $A\in\mathbb{R}^{I\times K}$4. This yields $A\in\mathbb{R}^{I\times K}$5 rather than $A\in\mathbb{R}^{I\times K}$6 for attaining $A\in\mathbb{R}^{I\times K}$7 (Khanna et al., 2017).

The centered self Khatri–Rao setting sharpens the picture for quadratic or covariance-type sensing models. If $A\in\mathbb{R}^{I\times K}$8 has iid columns drawn uniformly from the sphere of radius $A\in\mathbb{R}^{I\times K}$9 or with iid zero-mean, unit-variance sub-Gaussian entries, the self product has columns $B\in\mathbb{R}^{J\times K}$0. Because $B\in\mathbb{R}^{J\times K}$1, one defines the centered matrix $B\in\mathbb{R}^{J\times K}$2 with columns

$B\in\mathbb{R}^{J\times K}$3

and the normalized operator $B\in\mathbb{R}^{J\times K}$4 (Fengler et al., 2019). For this centered self Khatri–Rao product, there exist universal constants $B\in\mathbb{R}^{J\times K}$5 such that if

$B\in\mathbb{R}^{J\times K}$6

then

$B\in\mathbb{R}^{J\times K}$7

for any fixed $B\in\mathbb{R}^{J\times K}$8 (Fengler et al., 2019). 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 $B\in\mathbb{R}^{J\times K}$9 rather than in the vectors $A\odot B\in\mathbb{R}^{(IJ)\times K}$0 themselves (Fengler et al., 2019).

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 $A\odot B\in\mathbb{R}^{(IJ)\times K}$1 dense matrix (Sun et al., 2021). If $A\odot B\in\mathbb{R}^{(IJ)\times K}$2, choose $A\odot B\in\mathbb{R}^{(IJ)\times K}$3 and set

$A\odot B\in\mathbb{R}^{(IJ)\times K}$4

The memory cost is then $A\odot B\in\mathbb{R}^{(IJ)\times K}$5, compared with $A\odot B\in\mathbb{R}^{(IJ)\times K}$6 for a classical dense random projection and $A\odot B\in\mathbb{R}^{(IJ)\times K}$7 for a sparse projection with density $A\odot B\in\mathbb{R}^{(IJ)\times K}$8 (Sun et al., 2021). If $A\odot B\in\mathbb{R}^{(IJ)\times K}$9, this becomes $\ell$0, and for $\ell$1 with $\ell$2, storage is approximately $\ell$3 (Sun et al., 2021).

The basic isometry statement is unbiasedness: if each factor has independent columns in isotropic position, then

$\ell$4

for every $\ell$5 (Sun et al., 2021). When the entries of each factor are iid with mean $\ell$6, variance $\ell$7, and fourth moment $\ell$8, the variance of the squared embedding norm is

$\ell$9

A variance-reduced extension, denoted TRP(T), averages $A(:,\ell)\otimes B(:,\ell)$0 independent TRP maps and preserves the mean while reducing the first variance term to

$A(:,\ell)\otimes B(:,\ell)$1

(Sun et al., 2021).

For $A(:,\ell)\otimes B(:,\ell)$2, a non-asymptotic error bound is available under sub-Gaussian assumptions on the data vector $A(:,\ell)\otimes B(:,\ell)$3 and the factor entries: $A(:,\ell)\otimes B(:,\ell)$4 Via a standard JL-style argument, this yields an $A(:,\ell)\otimes B(:,\ell)$5-Johnson–Lindenstrauss embedding for $A(:,\ell)\otimes B(:,\ell)$6 points when

$A(:,\ell)\otimes B(:,\ell)$7

whereas dense Gaussian or fast JL methods require $A(:,\ell)\otimes B(:,\ell)$8 (Sun et al., 2021). The paper explicitly notes that the $A(:,\ell)\otimes B(:,\ell)$9 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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$00 memory instead of $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$01, and TRP(5) nearly matches the variance of dense random projections at approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$02 the storage for $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$03 (Sun et al., 2021). Sparse and very sparse choices of the factors preserve the “database-friendly” property of few nonzero queries and no floating-point operations (Sun et al., 2021).

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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$04. 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 (Saibaba et al., 31 Jul 2025, Bujanović et al., 2024). 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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$05, a random $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$06 is an $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$07-OSE if, with probability at least $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$08,

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$09

(Bujanović et al., 2024). For the two-factor Gaussian KRP

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$10

with iid $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$11 entries in the factors, one theorem states that $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$12 is an $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$13-OSE provided

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$14

(Bujanović et al., 2024). The dependence on $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$15 is therefore worse than the linear dependence of unstructured Gaussian sketches, but still polynomial rather than catastrophic.

A more general $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$16-factor subgaussian theory refines this picture. Let $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$17 have orthonormal columns and $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$18 be a KRP with i.i.d. mean-zero, unit-variance subgaussian entries in each factor. If

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$19

then with probability at least $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$20,

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$21

so $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$22 is a $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$23-subspace embedding for $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$24 (Saibaba et al., 31 Jul 2025). The same paper proves a Hanson–Wright-type tail bound for a single KRP column and a randomized range-finder guarantee: if $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$25 has rank-$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$26 leading singular space and the sketch dimension satisfies

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$27

then $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$28 obeys

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$29

with explicit $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$30 depending on $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$31, $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$32, $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$33, and $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$34 (Saibaba et al., 31 Jul 2025).

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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$35 and only poly-logarithmic dependence on the tensor order $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$36, rather than doubly exponential dependence (Saibaba et al., 31 Jul 2025). 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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$37 of the form

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$38

a single-view randomized SVD using KRPs draws

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$39

computes structured sketches

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$40

orthogonalizes $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$41, and recovers a core matrix via $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$42 (Saibaba et al., 31 Jul 2025). Assuming $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$43, the sketching cost is

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$44

with orthogonalization and least squares costing $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$45, compared with $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$46 for dense Gaussian sketches (Saibaba et al., 31 Jul 2025).

For Tucker compression, KRPs replace Gaussian sketches in RHOSVD and sequentially truncated HOSVD. In RHOSVD-KRP, each mode-$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$47 unfolding is multiplied by a KRP over all non-$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$48 modes,

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$49

followed by thin QR to obtain $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$50, after which the Tucker core is formed by multilinear projection (Saibaba et al., 31 Jul 2025). This naturally translates sketching into MTTKRP operations. Reported complexities are

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$51

for RHOSVD-KRP and

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$52

for RSTHOSVD-KRP under $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$53 and $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$54, with memoized RHOSVD-KRP further reducing MTTKRP cost by approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$55 (Saibaba et al., 31 Jul 2025).

In eigensolvers for Kronecker-structured operators, KRP sketches serve as structured starting blocks or projection matrices. For contour-integral methods, using $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$56 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 (Bujanović et al., 2024). In low-rank LOBPCG, if $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$57, matrix applications and preconditioner applications can be carried out in block-low-rank tensor format at cost $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$58 rather than $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$59 (Bujanović et al., 2024). The paper explicitly notes that LOBPCG is easier to integrate with KRP structure algorithmically, though with less theoretical justification than the contour-integral setting (Bujanović et al., 2024).

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 (Daas et al., 5 Nov 2025). In TT format, a $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$60-way tensor $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$61 is represented by cores

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$62

Rounding compresses these ranks, traditionally by deterministic TT-SVD or related procedures. The KRP-based approach instead sketches the large unfolding matrices using products

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$63

where each $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$64 is Gaussian (Daas et al., 5 Nov 2025).

The practical mechanism is a sequence of right-to-left partial contractions. For $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$65,

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$66

and for $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$67,

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$68

where $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$69 is the horizontal unfolding of the $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$70-th core (Daas et al., 5 Nov 2025). Each step is an MTTKRP costing $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$71, giving total cost

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$72

for the partial contractions and the same leading term overall for fixed-rank KRP-based TT-rounding (Daas et al., 5 Nov 2025).

The adaptive variant seeks $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$73 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

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$74

(Daas et al., 5 Nov 2025). The theoretical basis is a Frobenius-norm estimator: for

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$75

the quantity

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$76

is an unbiased estimator of $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$77, and if

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$78

then with probability at least $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$79,

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$80

(Daas et al., 5 Nov 2025). By summing the local errors in quadrature, the final TT approximation satisfies

$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$81

(Daas et al., 5 Nov 2025).

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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$82 over deterministic TT-rounding while maintaining comparable accuracy and only low overhead from adaptivity (Daas et al., 5 Nov 2025). 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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$83 of the random numbers of dense Gaussian sketching, runs approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$84 faster than dense Gaussian sketching and approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$85 faster than a block-Hankel-aware Gaussian method, at equal or better accuracy (Saibaba et al., 31 Jul 2025). In synthetic four-way Cauchy tensor compression, RHOSVD-KRP is approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$86–$B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$87 faster than RHOSVD at comparable approximation error, and the memoized version is the fastest among the tested methods (Saibaba et al., 31 Jul 2025). In tensor-based sensor placement for fluid flow, RHOSVD-KRP-MEMO is approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$88 faster than HOSVD and RHOSVD-KRP approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$89 faster, at essentially identical test error to RHOSVD (Saibaba et al., 31 Jul 2025). 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 $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$90 on $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$91 grids with per-iteration times reported between approximately $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$92 s and $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$93 s depending on the potential (Bujanović et al., 2024).

The literature also states several limitations. In TRP, variance scales with $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$94, so using too many factors degrades accuracy (Sun et al., 2021). The worst-case embedding dimension for TRP has an extra $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$95 factor in the $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$96 analysis, and the KRP OSE bound in eigensolver-oriented theory depends on $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$97 rather than linearly on $B=[b_1,\dots,b_n]\in\mathbb{R}^{m\times n}$98 (Sun et al., 2021, Bujanović et al., 2024). 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 (Bujanović et al., 2024). For TT rounding, the theoretical guarantees assume Gaussian factors and independence across cores (Daas et al., 5 Nov 2025).

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 (Khanna et al., 2017). 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 (Saibaba et al., 31 Jul 2025, Daas et al., 5 Nov 2025).

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 (Khanna et al., 2017, Sun et al., 2021, Bujanović et al., 2024, Saibaba et al., 31 Jul 2025, Daas et al., 5 Nov 2025).