Kronecker Product Least Squares

TL;DR

This paper introduces an iterative method using conjugate gradient and efficient Kronecker multiplication for solving Kronecker product least squares problems, demonstrating improved performance over direct methods.

Contribution

It presents a novel iterative approach combining conjugate gradient with Kronecker matrix-vector multiplication for least squares problems.

Findings

Iterative approach reduces computation time.

Method achieves comparable or better accuracy.

Numerical results validate efficiency improvements.

Abstract

In this rather brief note we present and discuss techniques for solving Kronecker matrix product least squares problems. Our main contribution is an iterative approach that uses the efficient Kronecker matrix-vector multiplication strategy (Fernandes et al. 1998) with a conjugate gradient solver. Numerical results contrast this approach---in terms of running times and accuracy---against the direct approach.

Tables

Table 1. Table 1: A comparison of the conjugate gradient solver and the full Kronecker product

No.	Size $({𝑨}_1,{𝑨}_2,{𝑨}_3,{𝑨}_4)$	Iterative method runtime (seconds)	Direct method runtime (seconds)
1	2,3,4,5	0.33	0.004
2	6,7,8,9	2.43	0.62
3	8,9,10,11	5.74	6.27
4	10,11,12,13	8.14	51.47

Equations

$$\underset{x}{\min }{\Vert \left({\mathbf{A}}_1\otimes {\mathbf{A}}_2\otimes \dots \otimes {\mathbf{A}}_n\right)\mathbf{x}-\mathbf{b}\Vert }_{l_2},$$

$$\mathbf{D}={\mathbf{A}}_1\otimes {\mathbf{A}}_2,$$

$$\left({\mathbf{A}}_1\otimes {\mathbf{A}}_2\right)\left({\mathbf{A}}_3\otimes {\mathbf{A}}_4\right)={\mathbf{A}}_1{\mathbf{A}}_3\otimes {\mathbf{A}}_2{\mathbf{A}}_4,$$

$$\left({\mathbf{A}}_1\otimes {\mathbf{A}}_2\right)\mathbf{x}={\mathbf{A}}_2\mathsf{mat}(\mathbf{x}){\mathbf{A}}_1^T$$

$${\left({\mathbf{A}}_1\otimes {\mathbf{A}}_2\right)}^T\left({\mathbf{A}}_1\otimes {\mathbf{A}}_2\right)\mathbf{x}={\left({\mathbf{A}}_1\otimes {\mathbf{A}}_2\right)}^T\mathbf{b}.$$

$$\left({\mathbf{R}}_1^T{\mathbf{Q}}_1^T\otimes {\mathbf{R}}_2^T{\mathbf{Q}}_2^T\right)\left({\mathbf{Q}}_1{\mathbf{R}}_1\otimes {\mathbf{Q}}_2{\mathbf{R}}_2\right)\mathbf{x}=\left({\mathbf{R}}_1^T{\mathbf{Q}}_1^T\otimes {\mathbf{R}}_2^T{\mathbf{Q}}_2^T\right)\mathbf{b}.$$

$$\left({\mathbf{R}}_1^T{\mathbf{R}}_1\otimes {\mathbf{R}}_2^T{\mathbf{R}}_2\right)\mathbf{x}=\left({\mathbf{R}}_1^T{\mathbf{Q}}_1^T\otimes {\mathbf{R}}_2^T{\mathbf{Q}}_2^T\right)\mathbf{b}$$

$${\mathbf{R}}_2^T{\mathbf{R}}_2\mathsf{mat}(\mathbf{x}){\mathbf{R}}_1^T{\mathbf{R}}_1={\mathbf{R}}_2^T{\mathbf{Q}}_2^T\mathsf{mat}(\mathbf{b}){\mathbf{Q}}_1{\mathbf{R}}_1,$$

$${\mathbf{R}}_2\mathsf{mat}(\mathbf{x}){\mathbf{R}}_1^T={\mathbf{Q}}_2^T\mathsf{mat}(\mathbf{b}){\mathbf{Q}}_1,$$

$$\Rightarrow \mathsf{mat}(\mathbf{x})={\mathbf{R}}_2^{-1}{\mathbf{Q}}_2^T\mathsf{mat}(\mathbf{b}){\mathbf{Q}}_1{\mathbf{R}}_1^{-T}.$$

$$\mathbf{x}\left({\otimes }_{k=1}^n{\mathbf{A}}_k\right)$$

$$=\mathbf{x}\prod _{k=1}^n{\mathbf{I}}_{(1:k-1)}\otimes {\mathbf{A}}_k\otimes {\mathbf{I}}_{(k+1:n)},$$

Full text

1 Introduction

Kronecker products routinely appear in uncertainty quantification [7], signal processing, image processing and more recently in quantum computing [8] and network science [6]. Their algebraic identities enable scalable and rapid matix computations, powering nuermous algorithms and solution techniques. Chief among these is the solution of linear systems, where Kronecker product identities have been applied to solving least squares [3], constrained least squares (using the null space method)[1], stacked least squares (using the generalized singular value decomposition; see page y of [5]) and weighted least squares [2]. Our objective in this manuscript is to offer solution strategies for the multilinear Kronecker product least squares problem

[TABLE]

where we assume that $\bm{A}_{k}$ for $k=1,\ldots,n$ are full rank square matrices. It is imperative that during the solution process to this linear system the full Kronecker product is never computed. While this problem has been addressed in literature for the case where $n=2$, to the best of our knowledge no solution exists for the more general case where $n>2$.

1.1 Notation and Preliminaries

We now introduce some notation and few Kronecker product properties that will be useful. Let $\bm{A}_{1}\in\mathbb{R}^{n\times n}$ and $\bm{A}_{2}\in\mathbb{R}^{m\times m}$, then their Kronecker product is

[TABLE]

where is a $\bm{D}\in\mathbb{R}^{mn\times mn}$ matrix. Essentially each block of $\bm{D}$ is formed by the product of $\bm{A}_{1}$ with the corresponding block of $\bm{A}_{2}$.

Kronecker products obey algebraic bilinearlity and associativity rules. We enumerate two key properties of Kronecker products that will be useful moving forwards.

The mixed product property: For four matrices $\bm{A}_{1},\bm{A}_{2},\bm{A}_{3}$ and $\bm{A}_{4}$

[TABLE]

assuming the sizes of the matrices permit the right hand side products to be formed. This property is especially useful in QR factorizations; 2. 2.

The Kronecker matrix vector product: For a vector $\mathbf{x}\in\mathbb{R}^{mn}$, the property

[TABLE]

holds. The operator $\mathsf{mat}(\mathbf{x})$ refers to a $\mathbb{R}^{m\times n}$ matrix comprising of all the $nm$ elements in the vector $\mathbf{x}$. Thus, one can interpret $\mathbf{x}$ as a vector of the stacked columns of $\mathsf{mat}(\mathbf{x})$.

These two rules are used to develop an expression for solving the Kronecker product least squares problem for the special case where $n=2$ in (1).

1.2 The Direct Least Squares Solution

Equation (1) may be re-written in the form of the normal equations

[TABLE]

Substituting matrices $\bm{A}_{1}$ and $\bm{A}_{2}$ with their thin QR factorizations gives us

[TABLE]

Now using the mixed-product property, (6) simplifies to

[TABLE]

[TABLE]

and since $\bm{R}_{1}$ and $\bm{R}_{2}$ are non-singular, this yields

[TABLE]

[TABLE]

Similar solution techniques may be achieved using a QR column pivoting factorization and the SVD (singular value decomposition) as discussed in [3]. Extending this approach to higher orders remains a challenge.

2 An Iterative Approach

In light of the above, we approach this problem using an iterative solution, which can in theory be generalized to any $n$. Our idea is to tie a conjugate gradient solver (see page 628 of [5]) with the effective matrix-vector product multiplication strategy in Fernandes et al.[4]; their strategy seeks to compute $\mathbf{x}\left(\otimes_{k=1}^{n}\bm{A}_{k}\right)$ without ever forming the full tensor product. Their algorithm is based on the tensor product decomposition

[TABLE]

where $\bm{I}$ are identity matrices. The notation $\bm{I}_{(1:k-1)}$ refers to an identity matrix with $\prod_{g=1}^{k-1}n_{g}$ rows and $\bm{I}_{(k+1:n)}$ to an identity matrix with $\prod_{g=k+1}^{n}n_{g}$ rows.

3 Numerical Results

We apply the iterative approach to the problem of approximating the coefficients of tensor grid polynomials. Our codes were implemented in MATLAB and utilized the software’s vanilla conjugate gradient solver. The experiments below were carried out on a 8GB Mac with a 2.8 GHz Intel Core i5 processor.

Table 1 presents representative results for the case where $n=4$ for different sizes for each of the matrices (which were all square). The table lists the total runtime for our iterative approach compared to that of the direct approach—which requires computing the full Kronecker product prior to solving the least squares problem.

In general it was observed that as the matrices get larger, the iterative technique proves to be more efficient from the perspective of run time. To evaluate the numerical accuracy of our approach, Figure 1 plots the least squares coefficients from a direct computation and from the iterative technique for case no. 4 in the table. Coefficient values above $10^{-8}$ (which was our tolerance limit for conjugate gradient solver), were found to match those of the direct solve.