∑MathIntermediate

Kronecker Product & Vec Operator

Key Points

•
The Kronecker product A ⊗ B expands a small matrix into a larger block matrix by multiplying every entry of A with the whole matrix B.
•
If A is m×n and B is p×q, then A ⊗ B is (mp)×(nq), arranged as blocks $a_{ij}$ B.
•
The vec operator stacks the columns of a matrix into one long column vector using column-major order.
•
A key identity is vec(AXB) = ( $B^{T}$ ⊗ A) vec(X), which lets you turn 2D matrix equations into 1D linear systems.
•
You should avoid explicitly forming A ⊗ B for large matrices; instead apply ( $B^{T}$ ⊗ A) to a vector via reshaping and usual multiplies.
•
Kronecker products preserve many properties: rank(A ⊗ B) = rank(A)rank(B), det(A ⊗ B) = det(A)^m det(B)^n, and eigenvalues multiply pairwise.
•
The Kronecker sum A ⊕ B = A ⊗ I + I ⊗ B models 2D grids and appears in PDEs and Sylvester equations.
•
In C++, careful indexing and consistent column-major vec are crucial to avoid transpose and ordering mistakes.

Prerequisites

→Matrix multiplication and transpose — Understanding standard matrix products and transposes is required to follow Kronecker definitions and the vec identity.
→Matrix indexing and storage order — Implementing vec/unvec and Kronecker correctly needs precise control of column-major stacking and indices.
→Solving linear systems — Vectorizing matrix equations leads to large linear systems that must be solved numerically.
→Eigenvalues and eigenvectors — Many Kronecker properties (spectrum of A ⊗ B and A ⊕ B) are phrased in spectral terms.
→Numerical stability and conditioning — Vectorized formulations can be ill-conditioned; understanding stability guides algorithm choice.

Detailed Explanation

Tap terms for definitions

01Overview

The Kronecker product is a way to build a very large, structured matrix from two smaller matrices. Given A and B, the Kronecker product A ⊗ B forms a block matrix whose blocks are just the entries of A scaling the entire matrix B. This operation shows up whenever multi-dimensional structure needs to be represented in a linear algebraic way: grids in numerical PDEs, separable filters in image processing, Markov chains on product spaces, and quantum computing with multi-qubit systems. The vec operator complements the Kronecker product. It converts a matrix into a long vector by stacking its columns on top of each other (column-major order). The best-known identity linking these two ideas is vec(AXB) = (B^T ⊗ A) vec(X), which converts a double-sided multiplication into a single large matrix acting on a vector. This identity is the workhorse behind many algorithms: it linearizes matrix equations like the Sylvester equation AX + XB = C and enables efficient application of Kronecker-structured operators without explicitly forming gigantic matrices. Intuitively, A ⊗ B captures how two independent linear transformations combine on a product space, while vec provides a bridge from 2D arrays to 1D vectors so we can apply standard linear algebra tools. Mastering these ideas yields simpler derivations, memory savings, and faster code for structured problems.

02Intuition & Analogies

Imagine you have two control knobs: one that scales and mixes rows (matrix A), and another that scales and mixes columns (matrix B). When you apply both knobs to a 2D grid (like a gray-scale image), the overall effect is separable: first adjust one direction, then the other. The Kronecker product is the big matrix that represents "doing both at once" on the flattened grid. Each entry a_{ij} of A becomes a whole block a_{ij}B because the column-direction transform B repeats for every row-direction interaction in A. That is why the size explodes: every row-column interaction in A duplicates the full behavior of B. Now picture the vec operator as folding a newspaper column by column: you take the first column top to bottom, then place the second column underneath it, and so on, until the entire page is one tall stack. Why do we do this? Because many algorithms, solvers, and libraries expect vectors. Vec lets us translate 2D matrix operations (like "left-multiply by A, right-multiply by B") into a single 1D operation using the Kronecker product. The golden rule is vec(AXB) = (B^T ⊗ A) vec(X): instead of computing AXB in 2D, you can flatten X and multiply by a big but structured matrix. Crucially, you rarely need to build that big matrix. If you reshape a vector back to a matrix X, then compute A X B with two standard multiplies, and re-stack with vec, you have effectively applied (B^T ⊗ A) to the vector without allocating the huge Kronecker matrix. This trick saves time and memory and is the key to practical use.

03Formal Definition

Let A ∈

R^{m \times n}

and B ∈

R^{p \times q}

. The Kronecker product A ⊗ B is the block matrix in

R^{(m p) \times (n q)}

defined by (A ⊗ B) = [

a_{ij}

B]_{

i = 1

..m,

j = 1

..n}, meaning its (i,j)-th block is the p×q matrix

a_{ij}

B. Index-wise, for 1 ≤

i \leq m

, 1 ≤

j \leq n

, 1 ≤

r \leq p

, 1 ≤

s \leq q

, (A ⊗ B)_{(i-1)p + r, (j-1)q + s} =

a_{ij}

b_{rs}

. The vec operator maps a matrix X ∈

R^{n \times p}

to a vector in

R^{n p}

by stacking columns: if X = [

x_{1}

x_{2}

...

x_{p}

] with

x_{k}

∈

R^{n}

, then vec(X) = [

x_{1}

;

x_{2}

; ...;

x_{p}

]. A central identity connecting the two is vec(AXB) = (

B^{⊤}

⊗ A) vec(X), which holds whenever the products are defined (A: m×n, X: n×p, B: p×q). Additional properties include the mixed-product rule (A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD) when dimensions match, and behavior under transpose and inverse: (A ⊗ B)^{

⊤

} =

A^{⊤}

⊗

B^{⊤}

, and if A and B are invertible, (A ⊗ B)^{-1} =

A^{- 1}

⊗

B^{- 1}

04When to Use

Use Kronecker products when modeling or exploiting separable structure across two independent dimensions.

Numerical PDEs and grids: Discrete 2D Laplacians on Cartesian grids are Kronecker sums A ⊕ B = A ⊗ I + I ⊗ B, enabling cleaner assembly and solvers.
Image processing and signals: Separable 2D filters (e.g., Gaussian blur) can be written as A X B, and applied efficiently via two 1D convolutions; vec identities help analyze and accelerate such operations.
Matrix equations: Linear equations in matrices like AX + XB = C can be vectorized to (I ⊗ A + B^{\top} ⊗ I) vec(X) = vec(C), allowing standard linear solvers.
Probabilistic models: Transition matrices of independent Markov chains combine via Kronecker products; stationary distributions and dynamics factor nicely.
Quantum computing and tensors: Multi-qubit gates and states are tensor products; in matrix form this is exactly the Kronecker product.
Large-scale linear algebra: To apply a huge operator with Kronecker structure to a vector, reshape the vector and use standard multiplies instead of forming the full matrix. In short, when your problem factors into two directions or components that act independently, Kronecker structure is often present and worth exploiting.

⚠️Common Mistakes

• Dimension blow-up: Explicitly forming A ⊗ B has size (mp)×(nq). Even modest inputs lead to huge outputs. Prefer applying (B^{\top} ⊗ A) to a vector by reshaping, computing A X B, and re-vectorizing. • Wrong vec ordering: This resource uses column-major vec. If you accidentally use row-major stacking, identities will need transposes or a commutation matrix. Be explicit about ordering in both code and math. • Missing transpose in the vec identity: The correct identity is vec(AXB) = (B^{\top} ⊗ A) vec(X). Forgetting the transpose on B is a common bug that flips or permutes results. • Confusing Kronecker with Hadamard: A ⊗ B is a block matrix (size multiplies). A ∘ B (Hadamard) is elementwise (size matches). Using the wrong one silently leads to dimension errors or nonsense results. • Inefficient multiplication: Multiplying by A ⊗ B via a generic dense mat-vec costs O(mnpq) and allocates huge memory. The reshape trick computes y = (B^{\top} ⊗ A) x as y = vec(A X B) in O(minimal) flops without building the big matrix. • Indexing off-by-ones: Implementations of kronecker and vec require careful index formulas. Unit tests on tiny matrices (2×2, 2×3) help catch mistakes early. • Numeric stability: Vectorizing matrix equations and solving the big system with naive Gaussian elimination can be less stable than specialized Sylvester solvers. Use appropriate algorithms for larger problems.

Key Formulas

Entry-wise Definition

(A \otimes B)_{(i - 1) p + r, (j - 1) q + s} = a_{ij} b_{rs}

Explanation: This formula defines each entry of the Kronecker product using block indexing. It says each $a_{ij}$ scales an entire copy of B placed in the appropriate block.

Block Definition

A \otimes B = a_{11} B ⋮ a_{m 1} B \dots ⋱ \dots a_{1 n} B ⋮ a_{mn} B

Explanation: The Kronecker product is a block matrix where the (i,j)-th block equals $a_{ij}$ times B. This visualizes how size multiplies with both dimensions.

Vec–Kronecker Identity

vec (A XB) = (B^{⊤} \otimes A) vec (X)

Explanation: This key identity turns a two-sided multiplication into a single large matrix times a vector. It enables solving matrix equations using standard linear solvers.

Mixed-Product Property

(A \otimes B) (C \otimes D) = (A C) \otimes (B D)

Explanation: Kronecker products interact nicely with normal products. When dimensions align, you can push products inside and multiply the smaller matrices.

Transpose Rule

(A \otimes B)^{⊤} = A^{⊤} \otimes B^{⊤}

Explanation: Transposing a Kronecker product equals the Kronecker product of transposes. This helps maintain symmetry and simplifies derivations.

Inverse Rule

(A \otimes B)^{- 1} = A^{- 1} \otimes B^{- 1}

Explanation: If A and B are invertible, so is their Kronecker product. The inverse distributes to each factor separately.

Rank Multiplicativity

rank (A \otimes B) = rank (A) rank (B)

Explanation: The rank of a Kronecker product is the product of the ranks. This is useful for reasoning about solvability and null spaces.

Determinant of Kronecker

det (A \otimes B) = det (A)^{m} det (B)^{n}

Explanation: For A ∈ $R^{n \times n}$ and B ∈ $R^{m \times m}$ , the determinant scales by powers of the other matrix's size. It describes volume scaling on product spaces.

Trace of Kronecker

tr (A \otimes B) = tr (A) tr (B)

Explanation: The trace factorizes across Kronecker products. It's handy in expectations and cost functions involving traces.

Eigenvalues of Kronecker

λ (A \otimes B) = {λ_{i} (A) μ_{j} (B)}_{i, j}

Explanation: Each eigenvalue of A ⊗ B is a product of an eigenvalue of A and an eigenvalue of B. Spectral properties separate across factors.

Frobenius Norm Multiplicativity

∥ A \otimes B ∥_{F} = ∥ A ∥_{F} ∥ B ∥_{F}

Explanation: The Frobenius norm of a Kronecker product equals the product of the Frobenius norms. Useful for error and complexity estimates.

Kronecker Sum

A \oplus B = A \otimes I_{m} + I_{n} \otimes B

Explanation: Defines the Kronecker sum for square matrices A and B. It models separable 2D operators such as the discrete Laplacian.

Eigenvalues of Kronecker Sum

λ (A \oplus B) = {λ_{i} (A) + μ_{j} (B)}_{i, j}

Explanation: The spectrum of the Kronecker sum consists of all pairwise sums of eigenvalues of A and B. This underpins fast solvers on grids.

Complexity Analysis

Let A be m×n and B be p×q. Forming the full Kronecker product

C = A

⊗ B requires writing an (mp)×(nq) matrix. Each of the m·n blocks is a scaled copy of B, which takes O(pq) time to fill, so overall time is O(mnpq). The space needed is also O(mnpq) to store all entries. This blow-up is the main practical limitation of explicit Kronecker construction. Applying a Kronecker-structured matrix to a vector can be done explicitly or implicitly. If you explicitly form K = (

B^{⊤}

⊗ A), then a matrix-vector multiply costs O(mnpq) time and O(mnpq) memory for K, which is prohibitive for large sizes. Instead, use the vec identity: reshape x ∈

R^{n p}

into X ∈

R^{n \times p}

, compute

Y = A

X B (order of operations chosen to minimize cost), and return vec(Y). If you compute

T = X

B in O(npq), followed by

Y = A

T in O(mnp), the total time is O(npq + mnp). Alternatively, compute

S = A

X in O(mnp), then

Y = S

B in O(mpq). Choose the cheaper path depending on relative sizes of m, n, p, q. The additional space beyond inputs and output is O(np + mq) for temporary matrices when done in-place-friendly fashion. For Kronecker sums

K = I_{p}

⊗ A +

B^{⊤}

⊗

I_{n}

used in vectorized Sylvester equations, forming K costs O(

n^{2}

p^{2}

) time/space for A, B ∈

R^{n \times n}

R^{p \times p}

. Solving

K v = c

via naive Gaussian elimination is O((np)^3) time and O((np)^2) space, which quickly becomes infeasible. In practice, specialized Sylvester solvers (Bartels–Stewart, Hessenberg–Schur) run in O(

n^{3}

p^{3}

) and avoid forming K, offering major speed and memory savings.

Code Examples

Kronecker product, vec/unvec, and verifying vec(AXB) = (B^T ⊗ A) vec(X)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Matrix {
5     size_t r, c;
6     vector<double> a; // row-major storage
7     Matrix(size_t r_=0, size_t c_=0, double val=0.0): r(r_), c(c_), a(r_*c_, val) {}
8     double& operator()(size_t i, size_t j) { return a[i*c + j]; }
9     double  operator()(size_t i, size_t j) const { return a[i*c + j]; }
10 };
11 
12 Matrix transpose(const Matrix& M){
13     Matrix T(M.c, M.r);
14     for(size_t i=0;i<M.r;++i)
15         for(size_t j=0;j<M.c;++j)
16             T(j,i) = M(i,j);
17     return T;
18 }
19 
20 Matrix multiply(const Matrix& A, const Matrix& B){
21     assert(A.c == B.r);
22     Matrix C(A.r, B.c, 0.0);
23     for(size_t i=0;i<A.r;++i){
24         for(size_t k=0;k<A.c;++k){
25             double aik = A(i,k);
26             for(size_t j=0;j<B.c;++j){
27                 C(i,j) += aik * B(k,j);
28             }
29         }
30     }
31     return C;
32 }
33 
34 Matrix kronecker(const Matrix& A, const Matrix& B){
35     Matrix K(A.r*B.r, A.c*B.c, 0.0);
36     for(size_t i=0;i<A.r;++i){
37         for(size_t j=0;j<A.c;++j){
38             double aij = A(i,j);
39             // place block aij * B at block position (i,j)
40             for(size_t r=0;r<B.r;++r){
41                 for(size_t s=0;s<B.c;++s){
42                     size_t R = i*B.r + r;
43                     size_t S = j*B.c + s;
44                     K(R,S) = aij * B(r,s);
45                 }
46             }
47         }
48     }
49     return K;
50 }
51 
52 // Column-major vec: stack columns of X
53 vector<double> vec_colmajor(const Matrix& X){
54     vector<double> v; v.reserve(X.r * X.c);
55     for(size_t j=0;j<X.c;++j)
56         for(size_t i=0;i<X.r;++i)
57             v.push_back(X(i,j));
58     return v;
59 }
60 
61 // Unvec (column-major): fill by columns to size rows x cols
62 Matrix unvec_colmajor(const vector<double>& v, size_t rows, size_t cols){
63     assert(v.size() == rows*cols);
64     Matrix X(rows, cols);
65     size_t t=0;
66     for(size_t j=0;j<cols;++j)
67         for(size_t i=0;i<rows;++i)
68             X(i,j) = v[t++];
69     return X;
70 }
71 
72 int main(){
73     // Small matrices to verify identity
74     Matrix A(2,2); A(0,0)=1; A(0,1)=2; A(1,0)=3; A(1,1)=4;      // 2x2
75     Matrix X(2,3); // 2x3
76     // Fill X with distinct values
77     int cnt=1; for(size_t i=0;i<X.r;++i) for(size_t j=0;j<X.c;++j) X(i,j)=cnt++;
78     Matrix B(3,2); // 3x2
79     // Fill B with small values
80     B(0,0)=1; B(0,1)=0; B(1,0)=-1; B(1,1)=2; B(2,0)=3; B(2,1)=-1;
81 
82     // Left side: vec(A X B)
83     Matrix AX = multiply(A, X);     // A*X  => (2x3)
84     Matrix AXB = multiply(AX, B);   // (2x2)
85     vector<double> lhs = vec_colmajor(AXB);
86 
87     // Right side: (B^T ⊗ A) vec(X)
88     Matrix BT = transpose(B);
89     Matrix K = kronecker(BT, A);    // ( (2x3) => 2x3 ) so K is (2*2)x(3*2) = 4x6
90     vector<double> vx = vec_colmajor(X); // length 2*3 = 6
91 
92     // Multiply K * vx
93     vector<double> rhs(K.r, 0.0);
94     for(size_t i=0;i<K.r;++i){
95         double sum=0;
96         for(size_t j=0;j<K.c;++j) sum += K(i,j) * vx[j];
97         rhs[i]=sum;
98     }
99 
100     cout.setf(ios::fixed); cout<<setprecision(3);
101     cout << "vec(A X B) = [";
102     for(double z: lhs) cout << z << ' '; cout << "]\n";
103     cout << "(B^T ⊗ A) vec(X) = [";
104     for(double z: rhs) cout << z << ' '; cout << "]\n";
105 
106     // Check maximum absolute difference
107     double err=0.0; for(size_t i=0;i<lhs.size();++i) err=max(err, fabs(lhs[i]-rhs[i]));
108     cout << "Max abs diff: " << err << "\n";
109     return 0;
110 }
111

This program implements a minimal dense Matrix type, Kronecker product, column-major vec/unvec, transpose, and standard multiplication. It verifies the core identity vec(AXB) = (B^T ⊗ A) vec(X) by computing both sides and comparing their entries.

Time: O(mnp + mpq + mnpq) for the demonstration (dominant term from forming Kronecker and doing K·vec).Space: O(mnpq) due to explicitly forming the Kronecker matrix.

Solving a tiny Sylvester equation AX + XB = C via Kronecker sum

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Matrix { size_t r,c; vector<double> a; Matrix(size_t r_=0,size_t c_=0,double v=0):r(r_),c(c_),a(r_*c_,v){} double& operator()(size_t i,size_t j){return a[i*c+j];} double operator()(size_t i,size_t j)const{return a[i*c+j];} };
5 
6 Matrix Iden(size_t n){ Matrix I(n,n,0.0); for(size_t i=0;i<n;++i) I(i,i)=1.0; return I; }
7 Matrix T(const Matrix& M){ Matrix R(M.c,M.r); for(size_t i=0;i<M.r;++i) for(size_t j=0;j<M.c;++j) R(j,i)=M(i,j); return R; }
8 Matrix mul(const Matrix& A,const Matrix& B){ assert(A.c==B.r); Matrix C(A.r,B.c,0.0); for(size_t i=0;i<A.r;++i) for(size_t k=0;k<A.c;++k){ double v=A(i,k); for(size_t j=0;j<B.c;++j) C(i,j)+=v*B(k,j);} return C; }
9 Matrix kron(const Matrix& A,const Matrix& B){ Matrix K(A.r*B.r,A.c*B.c,0.0); for(size_t i=0;i<A.r;++i) for(size_t j=0;j<A.c;++j){ double aij=A(i,j); for(size_t r=0;r<B.r;++r) for(size_t s=0;s<B.c;++s){ K(i*B.r+r, j*B.c+s) = aij*B(r,s); } } return K; }
10 vector<double> vec_col(const Matrix& X){ vector<double> v; v.reserve(X.r*X.c); for(size_t j=0;j<X.c;++j) for(size_t i=0;i<X.r;++i) v.push_back(X(i,j)); return v; }
11 Matrix unvec_col(const vector<double>& v,size_t r,size_t c){ assert(v.size()==r*c); Matrix X(r,c); size_t t=0; for(size_t j=0;j<c;++j) for(size_t i=0;i<r;++i) X(i,j)=v[t++]; return X; }
12 
13 // Simple Gaussian elimination with partial pivoting for dense systems
14 vector<double> solve_dense(Matrix A, vector<double> b){
15     size_t n=A.r; assert(A.r==A.c && b.size()==n);
16     // Augment
17     for(size_t i=0;i<n;++i) A.a.insert(A.a.begin()+ (i+1)*(A.c+i), b[i]); // not efficient; let's make explicit aug matrix
18     // We'll instead make explicit augmented matrix as 2D vector to keep it simple
19     vector<vector<double>> M(n, vector<double>(n+1));
20     for(size_t i=0;i<n;++i){ for(size_t j=0;j<n;++j) M[i][j]=A(i,j); M[i][n]=b[i]; }
21 
22     for(size_t col=0; col<n; ++col){
23         // pivot
24         size_t piv=col; for(size_t i=col+1;i<n;++i) if(fabs(M[i][col])>fabs(M[piv][col])) piv=i;
25         if(fabs(M[piv][col])<1e-12) throw runtime_error("Singular matrix");
26         swap(M[piv], M[col]);
27         double diag = M[col][col];
28         for(size_t j=col;j<=n;++j) M[col][j]/=diag; // normalize row
29         for(size_t i=0;i<n;++i){ if(i==col) continue; double f=M[i][col]; if(f==0) continue; for(size_t j=col;j<=n;++j) M[i][j]-=f*M[col][j]; }
30     }
31     vector<double> x(n); for(size_t i=0;i<n;++i) x[i]=M[i][n]; return x;
32 }
33 
34 int main(){
35     // A, B are 2x2; solve AX + X B = C
36     Matrix A(2,2); A(0,0)=2; A(0,1)=1; A(1,0)=-1; A(1,1)=3;
37     Matrix B(2,2); B(0,0)=1; B(0,1)=-2; B(1,0)=0; B(1,1)=4;
38     Matrix C(2,2); C(0,0)=1; C(0,1)=0; C(1,0)=2; C(1,1)=1;
39 
40     // Build K = I ⊗ A + B^T ⊗ I  (acting on vec(X))
41     Matrix I2 = Iden(2);
42     Matrix BT = T(B);
43     Matrix K = kron(I2, A); // I ⊗ A
44     Matrix BTkI = kron(BT, I2); // B^T ⊗ I
45     // K = K + BTkI
46     for(size_t i=0;i<K.r;++i) for(size_t j=0;j<K.c;++j) K(i,j) += BTkI(i,j);
47 
48     // Build right-hand side c = vec(C)
49     vector<double> cvec = vec_col(C);
50 
51     // Solve K v = c
52     vector<double> xvec = solve_dense(K, cvec);
53 
54     // Reshape into X and verify A X + X B ≈ C
55     Matrix X = unvec_col(xvec, 2, 2);
56     auto AX = mul(A, X);
57     auto XB = mul(X, B);
58     Matrix L(2,2); for(size_t i=0;i<2;++i) for(size_t j=0;j<2;++j) L(i,j)=AX(i,j)+XB(i,j);
59 
60     cout.setf(ios::fixed); cout<<setprecision(6);
61     cout << "Solution X:\n";
62     for(size_t i=0;i<2;++i){ for(size_t j=0;j<2;++j) cout << setw(12) << X(i,j); cout << '\n'; }
63     double err=0.0; for(size_t i=0;i<2;++i) for(size_t j=0;j<2;++j) err=max(err, fabs(L(i,j)-C(i,j)));
64     cout << "Max residual |AX+XB-C|_inf = " << err << '\n';
65     return 0;
66 }
67

This example vectorizes the Sylvester equation AX + XB = C into a linear system (I ⊗ A + B^T ⊗ I) vec(X) = vec(C), builds the Kronecker-sum matrix, solves the tiny dense system with Gaussian elimination, and reshapes the solution back to X. It then verifies the residual.

Time: O(n^2 p^2) to build K plus O((np)^3) to solve with naive elimination (here n=p=2).Space: O(n^2 p^2) for K and O((np)^2) during elimination.

Apply y = (B^T ⊗ A) x without forming the Kronecker matrix

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Matrix{ size_t r,c; vector<double>a; Matrix(size_t r_=0,size_t c_=0,double v=0):r(r_),c(c_),a(r_*c_,v){} double& operator()(size_t i,size_t j){return a[i*c+j];} double operator()(size_t i,size_t j)const{return a[i*c+j];} };
5 Matrix mul(const Matrix& A,const Matrix& B){ assert(A.c==B.r); Matrix C(A.r,B.c,0.0); for(size_t i=0;i<A.r;++i) for(size_t k=0;k<A.c;++k){ double v=A(i,k); for(size_t j=0;j<B.c;++j) C(i,j)+=v*B(k,j);} return C; }
6 Matrix T(const Matrix& M){ Matrix R(M.c,M.r); for(size_t i=0;i<M.r;++i) for(size_t j=0;j<M.c;++j) R(j,i)=M(i,j); return R; }
7 vector<double> vec_col(const Matrix& X){ vector<double> v; v.reserve(X.r*X.c); for(size_t j=0;j<X.c;++j) for(size_t i=0;i<X.r;++i) v.push_back(X(i,j)); return v; }
8 Matrix unvec_col(const vector<double>& v,size_t r,size_t c){ assert(v.size()==r*c); Matrix X(r,c); size_t t=0; for(size_t j=0;j<c;++j) for(size_t i=0;i<r;++i) X(i,j)=v[t++]; return X; }
9 
10 // Efficient apply: y = (B^T ⊗ A) x by reshaping x to X (n×p), computing Y = A X B, returning vec(Y)
11 vector<double> apply_kron_matvec(const Matrix& A, const Matrix& B, const vector<double>& x, size_t n, size_t p){
12     // A: m×n, B: p×q, x length n*p, interpret x as vec(X) with X ∈ R^{n×p}
13     size_t m = A.r, q = B.c;
14     assert(A.c==n && B.r==p && x.size()==n*p);
15     Matrix X = unvec_col(x, n, p);
16     Matrix XB = mul(X, B);    // n×q
17     Matrix Y  = mul(A, XB);   // m×q
18     return vec_col(Y);        // length m*q
19 }
20 
21 int main(){
22     // Small check against explicit construction for correctness
23     size_t m=3,n=2,p=2,q=3;
24     Matrix A(m,n), B(p,q);
25     // Fill A,B with simple values
26     for(size_t i=0;i<m;++i) for(size_t j=0;j<n;++j) A(i,j)= (i+1) + 0.1*(j+1);
27     for(size_t i=0;i<p;++i) for(size_t j=0;j<q;++j) B(i,j)= (j+1) - 0.2*(i+1);
28 
29     // Random X and x=vec(X)
30     Matrix X(n,p); std::mt19937 rng(42); std::uniform_real_distribution<double> U(-1,1);
31     for(size_t i=0;i<n;++i) for(size_t j=0;j<p;++j) X(i,j)=U(rng);
32     vector<double> xvec = vec_col(X);
33 
34     // Fast apply
35     vector<double> y_fast = apply_kron_matvec(A, B, xvec, n, p);
36 
37     // Explicit check via forming K = (B^T ⊗ A)
38     // Build Kronecker explicitly (for testing only)
39     auto kron = [&](const Matrix& M,const Matrix& N){ Matrix K(M.r*N.r, M.c*N.c); for(size_t i=0;i<M.r;++i) for(size_t j=0;j<M.c;++j){ double mij=M(i,j); for(size_t r=0;r<N.r;++r) for(size_t s=0;s<N.c;++s) K(i*N.r+r, j*N.c+s)=mij*N(r,s);} return K; };
40     Matrix BT = T(B);
41     Matrix K = kron(BT, A);
42     vector<double> y_exp(K.r, 0.0);
43     for(size_t i=0;i<K.r;++i){ double s=0; for(size_t j=0;j<K.c;++j) s+=K(i,j)*xvec[j]; y_exp[i]=s; }
44 
45     // Compare
46     double err=0.0; for(size_t i=0;i<y_exp.size();++i) err=max(err, fabs(y_exp[i]-y_fast[i]));
47     cout.setf(ios::fixed); cout<<setprecision(9);
48     cout << "Max abs diff (fast vs explicit): " << err << "\n";
49     cout << "Output length: " << y_fast.size() << " (should be m*q = " << m*q << ")\n";
50     return 0;
51 }
52

This code applies y = (B^T ⊗ A) x using the reshape trick: interpret x as vec(X), compute Y = A X B, and return vec(Y). It verifies correctness against an explicit Kronecker construction on small sizes. In real applications, you should always prefer the fast method to avoid the O(mnpq) memory blow-up.

Time: O(n p q + m n q) using order (X B) then A(XB), or O(m n p + m p q) with the alternative order.Space: O(n p + m q) extra for temporaries, versus O(m n p q) if forming the Kronecker matrix.

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct Matrix {
5	size_t r, c;
6	vector<double> a; // row-major storage
7	Matrix(size_t r_=0, size_t c_=0, double val=0.0): r(r_), c(c_), a(r_*c_, val) {}
8	double& operator()(size_t i, size_t j) { return a[i*c + j]; }
9	double operator()(size_t i, size_t j) const { return a[i*c + j]; }
10	};
11
12	Matrix transpose(const Matrix& M){
13	Matrix T(M.c, M.r);
14	for(size_t i=0;i<M.r;++i)
15	for(size_t j=0;j<M.c;++j)
16	T(j,i) = M(i,j);
17	return T;
18	}
19
20	Matrix multiply(const Matrix& A, const Matrix& B){
21	assert(A.c == B.r);
22	Matrix C(A.r, B.c, 0.0);
23	for(size_t i=0;i<A.r;++i){
24	for(size_t k=0;k<A.c;++k){
25	double aik = A(i,k);
26	for(size_t j=0;j<B.c;++j){
27	C(i,j) += aik * B(k,j);
28	}
29	}
30	}
31	return C;
32	}
33
34	Matrix kronecker(const Matrix& A, const Matrix& B){
35	Matrix K(A.rB.r, A.cB.c, 0.0);
36	for(size_t i=0;i<A.r;++i){
37	for(size_t j=0;j<A.c;++j){
38	double aij = A(i,j);
39	// place block aij * B at block position (i,j)
40	for(size_t r=0;r<B.r;++r){
41	for(size_t s=0;s<B.c;++s){
42	size_t R = i*B.r + r;
43	size_t S = j*B.c + s;
44	K(R,S) = aij * B(r,s);
45	}
46	}
47	}
48	}
49	return K;
50	}
51
52	// Column-major vec: stack columns of X
53	vector<double> vec_colmajor(const Matrix& X){
54	vector<double> v; v.reserve(X.r * X.c);
55	for(size_t j=0;j<X.c;++j)
56	for(size_t i=0;i<X.r;++i)
57	v.push_back(X(i,j));
58	return v;
59	}
60
61	// Unvec (column-major): fill by columns to size rows x cols
62	Matrix unvec_colmajor(const vector<double>& v, size_t rows, size_t cols){
63	assert(v.size() == rows*cols);
64	Matrix X(rows, cols);
65	size_t t=0;
66	for(size_t j=0;j<cols;++j)
67	for(size_t i=0;i<rows;++i)
68	X(i,j) = v[t++];
69	return X;
70	}
71
72	int main(){
73	// Small matrices to verify identity
74	Matrix A(2,2); A(0,0)=1; A(0,1)=2; A(1,0)=3; A(1,1)=4; // 2x2
75	Matrix X(2,3); // 2x3
76	// Fill X with distinct values
77	int cnt=1; for(size_t i=0;i<X.r;++i) for(size_t j=0;j<X.c;++j) X(i,j)=cnt++;
78	Matrix B(3,2); // 3x2
79	// Fill B with small values
80	B(0,0)=1; B(0,1)=0; B(1,0)=-1; B(1,1)=2; B(2,0)=3; B(2,1)=-1;
81
82	// Left side: vec(A X B)
83	Matrix AX = multiply(A, X); // A*X => (2x3)
84	Matrix AXB = multiply(AX, B); // (2x2)
85	vector<double> lhs = vec_colmajor(AXB);
86
87	// Right side: (B^T ⊗ A) vec(X)
88	Matrix BT = transpose(B);
89	Matrix K = kronecker(BT, A); // ( (2x3) => 2x3 ) so K is (22)x(32) = 4x6
90	vector<double> vx = vec_colmajor(X); // length 2*3 = 6
91
92	// Multiply K * vx
93	vector<double> rhs(K.r, 0.0);
94	for(size_t i=0;i<K.r;++i){
95	double sum=0;
96	for(size_t j=0;j<K.c;++j) sum += K(i,j) * vx[j];
97	rhs[i]=sum;
98	}
99
100	cout.setf(ios::fixed); cout<<setprecision(3);
101	cout << "vec(A X B) = [";
102	for(double z: lhs) cout << z << ' '; cout << "]\n";
103	cout << "(B^T ⊗ A) vec(X) = [";
104	for(double z: rhs) cout << z << ' '; cout << "]\n";
105
106	// Check maximum absolute difference
107	double err=0.0; for(size_t i=0;i<lhs.size();++i) err=max(err, fabs(lhs[i]-rhs[i]));
108	cout << "Max abs diff: " << err << "\n";
109	return 0;
110	}
111

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct Matrix { size_t r,c; vector<double> a; Matrix(size_t r_=0,size_t c_=0,double v=0):r(r_),c(c_),a(r_c_,v){} double& operator()(size_t i,size_t j){return a[ic+j];} double operator()(size_t i,size_t j)const{return a[i*c+j];} };
5
6	Matrix Iden(size_t n){ Matrix I(n,n,0.0); for(size_t i=0;i<n;++i) I(i,i)=1.0; return I; }
7	Matrix T(const Matrix& M){ Matrix R(M.c,M.r); for(size_t i=0;i<M.r;++i) for(size_t j=0;j<M.c;++j) R(j,i)=M(i,j); return R; }
8	Matrix mul(const Matrix& A,const Matrix& B){ assert(A.c==B.r); Matrix C(A.r,B.c,0.0); for(size_t i=0;i<A.r;++i) for(size_t k=0;k<A.c;++k){ double v=A(i,k); for(size_t j=0;j<B.c;++j) C(i,j)+=v*B(k,j);} return C; }
9	Matrix kron(const Matrix& A,const Matrix& B){ Matrix K(A.rB.r,A.cB.c,0.0); for(size_t i=0;i<A.r;++i) for(size_t j=0;j<A.c;++j){ double aij=A(i,j); for(size_t r=0;r<B.r;++r) for(size_t s=0;s<B.c;++s){ K(iB.r+r, jB.c+s) = aij*B(r,s); } } return K; }
10	vector<double> vec_col(const Matrix& X){ vector<double> v; v.reserve(X.r*X.c); for(size_t j=0;j<X.c;++j) for(size_t i=0;i<X.r;++i) v.push_back(X(i,j)); return v; }
11	Matrix unvec_col(const vector<double>& v,size_t r,size_t c){ assert(v.size()==r*c); Matrix X(r,c); size_t t=0; for(size_t j=0;j<c;++j) for(size_t i=0;i<r;++i) X(i,j)=v[t++]; return X; }
12
13	// Simple Gaussian elimination with partial pivoting for dense systems
14	vector<double> solve_dense(Matrix A, vector<double> b){
15	size_t n=A.r; assert(A.r==A.c && b.size()==n);
16	// Augment
17	for(size_t i=0;i<n;++i) A.a.insert(A.a.begin()+ (i+1)*(A.c+i), b[i]); // not efficient; let's make explicit aug matrix
18	// We'll instead make explicit augmented matrix as 2D vector to keep it simple
19	vector<vector<double>> M(n, vector<double>(n+1));
20	for(size_t i=0;i<n;++i){ for(size_t j=0;j<n;++j) M[i][j]=A(i,j); M[i][n]=b[i]; }
21
22	for(size_t col=0; col<n; ++col){
23	// pivot
24	size_t piv=col; for(size_t i=col+1;i<n;++i) if(fabs(M[i][col])>fabs(M[piv][col])) piv=i;
25	if(fabs(M[piv][col])<1e-12) throw runtime_error("Singular matrix");
26	swap(M[piv], M[col]);
27	double diag = M[col][col];
28	for(size_t j=col;j<=n;++j) M[col][j]/=diag; // normalize row
29	for(size_t i=0;i<n;++i){ if(i==col) continue; double f=M[i][col]; if(f==0) continue; for(size_t j=col;j<=n;++j) M[i][j]-=f*M[col][j]; }
30	}
31	vector<double> x(n); for(size_t i=0;i<n;++i) x[i]=M[i][n]; return x;
32	}
33
34	int main(){
35	// A, B are 2x2; solve AX + X B = C
36	Matrix A(2,2); A(0,0)=2; A(0,1)=1; A(1,0)=-1; A(1,1)=3;
37	Matrix B(2,2); B(0,0)=1; B(0,1)=-2; B(1,0)=0; B(1,1)=4;
38	Matrix C(2,2); C(0,0)=1; C(0,1)=0; C(1,0)=2; C(1,1)=1;
39
40	// Build K = I ⊗ A + B^T ⊗ I (acting on vec(X))
41	Matrix I2 = Iden(2);
42	Matrix BT = T(B);
43	Matrix K = kron(I2, A); // I ⊗ A
44	Matrix BTkI = kron(BT, I2); // B^T ⊗ I
45	// K = K + BTkI
46	for(size_t i=0;i<K.r;++i) for(size_t j=0;j<K.c;++j) K(i,j) += BTkI(i,j);
47
48	// Build right-hand side c = vec(C)
49	vector<double> cvec = vec_col(C);
50
51	// Solve K v = c
52	vector<double> xvec = solve_dense(K, cvec);
53
54	// Reshape into X and verify A X + X B ≈ C
55	Matrix X = unvec_col(xvec, 2, 2);
56	auto AX = mul(A, X);
57	auto XB = mul(X, B);
58	Matrix L(2,2); for(size_t i=0;i<2;++i) for(size_t j=0;j<2;++j) L(i,j)=AX(i,j)+XB(i,j);
59
60	cout.setf(ios::fixed); cout<<setprecision(6);
61	cout << "Solution X:\n";
62	for(size_t i=0;i<2;++i){ for(size_t j=0;j<2;++j) cout << setw(12) << X(i,j); cout << '\n'; }
63	double err=0.0; for(size_t i=0;i<2;++i) for(size_t j=0;j<2;++j) err=max(err, fabs(L(i,j)-C(i,j)));
64	cout << "Max residual \|AX+XB-C\|_inf = " << err << '\n';
65	return 0;
66	}
67

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct Matrix{ size_t r,c; vector<double>a; Matrix(size_t r_=0,size_t c_=0,double v=0):r(r_),c(c_),a(r_c_,v){} double& operator()(size_t i,size_t j){return a[ic+j];} double operator()(size_t i,size_t j)const{return a[i*c+j];} };
5	Matrix mul(const Matrix& A,const Matrix& B){ assert(A.c==B.r); Matrix C(A.r,B.c,0.0); for(size_t i=0;i<A.r;++i) for(size_t k=0;k<A.c;++k){ double v=A(i,k); for(size_t j=0;j<B.c;++j) C(i,j)+=v*B(k,j);} return C; }
6	Matrix T(const Matrix& M){ Matrix R(M.c,M.r); for(size_t i=0;i<M.r;++i) for(size_t j=0;j<M.c;++j) R(j,i)=M(i,j); return R; }
7	vector<double> vec_col(const Matrix& X){ vector<double> v; v.reserve(X.r*X.c); for(size_t j=0;j<X.c;++j) for(size_t i=0;i<X.r;++i) v.push_back(X(i,j)); return v; }
8	Matrix unvec_col(const vector<double>& v,size_t r,size_t c){ assert(v.size()==r*c); Matrix X(r,c); size_t t=0; for(size_t j=0;j<c;++j) for(size_t i=0;i<r;++i) X(i,j)=v[t++]; return X; }
9
10	// Efficient apply: y = (B^T ⊗ A) x by reshaping x to X (n×p), computing Y = A X B, returning vec(Y)
11	vector<double> apply_kron_matvec(const Matrix& A, const Matrix& B, const vector<double>& x, size_t n, size_t p){
12	// A: m×n, B: p×q, x length n*p, interpret x as vec(X) with X ∈ R^{n×p}
13	size_t m = A.r, q = B.c;
14	assert(A.c==n && B.r==p && x.size()==n*p);
15	Matrix X = unvec_col(x, n, p);
16	Matrix XB = mul(X, B); // n×q
17	Matrix Y = mul(A, XB); // m×q
18	return vec_col(Y); // length m*q
19	}
20
21	int main(){
22	// Small check against explicit construction for correctness
23	size_t m=3,n=2,p=2,q=3;
24	Matrix A(m,n), B(p,q);
25	// Fill A,B with simple values
26	for(size_t i=0;i<m;++i) for(size_t j=0;j<n;++j) A(i,j)= (i+1) + 0.1*(j+1);
27	for(size_t i=0;i<p;++i) for(size_t j=0;j<q;++j) B(i,j)= (j+1) - 0.2*(i+1);
28
29	// Random X and x=vec(X)
30	Matrix X(n,p); std::mt19937 rng(42); std::uniform_real_distribution<double> U(-1,1);
31	for(size_t i=0;i<n;++i) for(size_t j=0;j<p;++j) X(i,j)=U(rng);
32	vector<double> xvec = vec_col(X);
33
34	// Fast apply
35	vector<double> y_fast = apply_kron_matvec(A, B, xvec, n, p);
36
37	// Explicit check via forming K = (B^T ⊗ A)
38	// Build Kronecker explicitly (for testing only)
39	auto kron = [&](const Matrix& M,const Matrix& N){ Matrix K(M.rN.r, M.cN.c); for(size_t i=0;i<M.r;++i) for(size_t j=0;j<M.c;++j){ double mij=M(i,j); for(size_t r=0;r<N.r;++r) for(size_t s=0;s<N.c;++s) K(iN.r+r, jN.c+s)=mij*N(r,s);} return K; };
40	Matrix BT = T(B);
41	Matrix K = kron(BT, A);
42	vector<double> y_exp(K.r, 0.0);
43	for(size_t i=0;i<K.r;++i){ double s=0; for(size_t j=0;j<K.c;++j) s+=K(i,j)*xvec[j]; y_exp[i]=s; }
44
45	// Compare
46	double err=0.0; for(size_t i=0;i<y_exp.size();++i) err=max(err, fabs(y_exp[i]-y_fast[i]));
47	cout.setf(ios::fixed); cout<<setprecision(9);
48	cout << "Max abs diff (fast vs explicit): " << err << "\n";
49	cout << "Output length: " << y_fast.size() << " (should be mq = " << mq << ")\n";
50	return 0;
51	}
52