# Solving linear system with sum of kronecker products

**URL:** <https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253>\
**Category:** Numerics\
**Tags:** question, linearalgebra\
**Created:** [January 25, 2024, 1:11pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253 "2024-01-25T13:11:52Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 25, 2024, 1:11pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/1 "2024-01-25T13:11:53Z")

</div>

Suppose I have

((A \otimes B) + (C \otimes D))x = v

where A and C and B and D have the same size, respectively. Suppose also that I particularly want to avoid instantiating the Kronecker products since they are large.

What would be the recommended way to solve for `x`? My default would be using IterativeSolvers.jl, using Kronecker.jl to represent the \otimes products. But I am not sure what to use for the sum, other than `LazyArrays.BroadcastArray`.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 25, 2024, 1:59pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/2 "2024-01-25T13:59:56Z")

</div>

Assuming A, B, C, and D are square, which seems reasonable since you could otherwise have an over or under determined system, the usual thing to do would be to use Kronecker product properties to rewrite this as a matrix equation

BXA^T+ DXC^T = V

where x=\mbox{vec}(X) and v = \mbox{vec}(V). Here \mbox{vec} stacks the columns of a matrix into a vector. This is a generalized Sylvester equation. There are then solvers that compute generalized Schur decompositions of (B,D) and (A,C) to allow for the efficient direct column-wise solution of X using a back-substitution procedure. If every matrix is n\times n, it’s an O(n^3) algorithm to compute the n^2 elements of X. I believe that [MatrixEquations.jl](https://github.com/andreasvarga/MatrixEquations.jl) has this implemented with the function `gsylv`.

Even if your matrices are sparse there is limited benefit in trying to use an iterative method to exploit sparsity, since X will not in general be sparse. The direct method should be the way to go.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [January 25, 2024, 3:02pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/3 "2024-01-25T15:02:03Z")

</div>

To complete the answer, if you still wanna use an iterative solver, you don’t need to “represent” the sum in any way other than a functional form

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 25, 2024, 3:20pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/4 "2024-01-25T15:20:07Z")

</div>

Thanks.

Some background: I am solving a Fokker-Planck PDE using collocation. Suppose that T\_0, T\_1 are the collocation matrices for values and derivatives of a basis for time, and S\_0, S\_1, S\_2 are similarly for a state, and c\_0, etc are various coefficients.

Then generally my left hand side matrix is something like

T\_1 ⊗ S\_0 + T\_0 ⊗ (\mathrm{diag}(c\_1) S\_0) + T\_0 ⊗ (\mathrm{diag}(c\_0) S\_1) + \dots

which are for the drift only, T\_0 \otimes S\_2 etc terms follow if I have diffusion.

Then

D = \mathrm{diag}(c\_1) S\_0 + \mathrm{diag}(c\_0) S\_1 + \dots

which makes my life very, very simple and my computation fast.

Even if the matrix is not sparse, the fact that I don’t have to form the Kronecker explicitly is a huge time saver.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 25, 2024, 4:27pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/5 "2024-01-25T16:27:22Z")

</div>

I did a double take seeing the sum of more than two Kronecker products before I saw your point about T\_0 being common to all of them except the first. If you have 3 or more Kronecker products, you don’t have the nice direct method and, as far as I know, you have to use iterative methods. As you say, your problem is a lot nicer with the way you are able to form D.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 25, 2024, 4:31pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/6 "2024-01-25T16:31:27Z")

</div>

> But I am not sure what to use for the sum, other than `LazyArrays.BroadcastArray` .

I am not sure I understand your problem. Can you evaluate (A\otimes B) x and (C\otimes D) x ?

Then, looking at the above comment regarding the FP equation, the convergence could be slow and preconditioning by the Laplacian would help i.m.o.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 25, 2024, 4:45pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/7 "2024-01-25T16:45:39Z")

</div>

> [@rveltz](#):
>
> Can you evaluate

yes

> [@rveltz](#):
>
> preconditioning by the Laplacian would help i.m.o.

Thanks. I am new to PDE solving and not familiar with the concept, can you please recommend a text that discusses preconditioning?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 25, 2024, 5:16pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/8 "2024-01-25T17:16:54Z")

</div>

Basically, solving

$$\Delta u + u = v$$

with gmres has poor convergence. However, solving

$$u + \Delta^{-1}u = \Delta^{-1}v$$

is fast. You can pass the preconditioner `Pl=lu(Delta)` to gmres and it will do the above conversion for you. See [Restarted GMRES · IterativeSolvers.jl](https://iterativesolvers.julialinearalgebra.org/stable/linear_systems/gmres/)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 25, 2024, 5:19pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/9 "2024-01-25T17:19:06Z")

</div>

> [@Tamas\_Papp](#):
>
> yes

Then you can pass the function written in a lazy way

```julia
function LinOp(x)
    res = kron(A,B)*x + kron(C,D)*x
end

```

to `gmres`

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 26, 2024, 10:08am UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/10 "2024-01-26T10:08:19Z")

</div>

Thanks. I found that my problem is reasonably well-conditioned because the spectral bases are (rationally transformed) Chebyshev polynomials, so I did not look for a preconditioner.

Is the solution from `gmres` supposed to be close to `A \ b`, or am I using it in the wrong way? Because for my problem, it is wildly off, and so is the residual:

```julia
julia> using IterativeSolvers

julia> θ1 = A \ b;

julia> θ2 = gmres(A, b);

julia> maximum(abs, A * θ1 - b)
5.551115123125783e-16

julia> maximum(abs, A * θ2 - b)
0.5582449534355522

```

It this is interesting to the devs I can `collect` `A` and `b` as a matrix (as I did above) and submit an MWE as an issue.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 26, 2024, 12:09pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/11 "2024-01-26T12:09:49Z")

</div>

you can check for convergence by doing

```julia
gmres(A, b; log = true, verbose = true)

```

> Thanks. I found that my problem is reasonably well-conditioned because the spectral bases are (rationally transformed) Chebyshev polynomials, so I did not look for a preconditioner.

That contradicts the outcome of your numerical experiment it seems

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 26, 2024, 12:52pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/12 "2024-01-26T12:52:35Z")

</div>

> [@Tamas\_Papp](#):
>
> Thanks. I found that my problem is reasonably well-conditioned because the spectral bases are (rationally transformed) Chebyshev polynomials, so I did not look for a preconditioner.

I’m probably going out on a limb since I have no experience with Fokker-Planck equations, but it might be worth noting that convergence with GMRES defies simple characterization. It can perform badly even with a well conditioned matrix.  
Preconditioning might still be necessary and, unless there’s some bug here, it looks like it probably is necessary for your problem.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [January 26, 2024, 12:57pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/13 "2024-01-26T12:57:13Z")

</div>

[A\_X \otimes B\_Y + C\_X \otimes D\_Y]x=v  
=\>A\_X^{-1} \otimes D\_Y^{-1}[A\_X \otimes B\_Y + C\_X \otimes D\_Y]x=A\_X^{-1} \otimes D\_Y^{-1}v  
=\> [1\_X \otimes D\_Y^{-1}B\_Y + A\_X^{-1} C\_X \otimes 1\_Y]x=A\_X^{-1} \otimes D\_Y^{-1}v  
let w=A\_X^{-1} \otimes D\_Y^{-1}v , E\_Y=D\_Y^{-1}B\_Y and E\_X=A\_X^{-1} C\_X  
=\>[E\_X \otimes 1\_Y + 1\_X \otimes E\_Y]x=w  
that can be effciently solved with a direct method via the “tensor trick” :

Let E\_X=M\_X \Pi\_XM\_X^{-1} and E\_Y=M\_Y \Pi\_YM\_Y^{-1} the diagonalization of E\_X and E\_Y.

=\>[E\_X \otimes 1\_Y + 1\_X \otimes E\_Y]x=w  
P\_{XY} = M\_X\otimes 1\_Y + 1\_X\otimes M\_Y is diagonal and can be simply inverted  
[P^{-1}\_{XY}]\_{\alpha\beta,\alpha'\beta'} = \frac{\delta\_{\alpha,\alpha'},\delta\_{\beta,\beta'}}{l\_{x\alpha}+l\_{y\beta}} where l\_{x\alpha} and l\_{y\beta} are the \Pi\_X and \Pi\_Y diagonal elements.  
[E\_X \otimes 1\_Y + 1\_X \otimes E\_Y]x=[M\_X \otimes M\_Y]P\_{XY}[M^{-1}\_X \otimes M\_Y^{-1}]x=w  
=\>x=[M\_X \otimes M\_Y]P^{-1}\_{XY}[M^{-1}\_X \otimes M\_Y^{-1}]w a direct computation involving only small matrices (size of A\_X).

I did implement the tensor trick with Julia in `src/poisson2D_TT.jl` here [GitHub - triscale-innov/LidJul.jl: GMG, Poisson solver and Lid cavity](https://github.com/triscale-innov/LidJul.jl)

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 26, 2024, 1:19pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/14 "2024-01-26T13:19:22Z")

</div>

> [@mstewart](#):
>
> Preconditioning might still be necessary

Sorry, I still only have a vague clue about preconditioning despite reading up about it in the context of some methods (particularly CG).

Is there a more robust method that is “matrix-free”, in the sense the evaluating Ax is sufficient?

> [@LaurentPlagne](#):
>
> that can be effciently solved with a direct method via the “tensor trick”

Can you please explain or provide a link (does a package do this? Kronecker.jl seems to provide a type like this but no `\` method as far as I can tell).

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [January 26, 2024, 1:51pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/15 "2024-01-26T13:51:21Z")

</div>

I did edit my previous post. You can have a look [https://www.sciencedirect.com/science/article/abs/pii/S0021999199963386](https://www.sciencedirect.com/science/article/abs/pii/S0021999199963386)

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 26, 2024, 2:08pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/16 "2024-01-26T14:08:08Z")

</div>

> [@Tamas\_Papp](#):
>
> Sorry, I still only have a vague clue about preconditioning despite reading up about it in the context of some methods (particularly CG).
> 
> Is there a more robust method that is “matrix-free”, in the sense the evaluating Ax is sufficient?

Maybe someone who is more of an iterative methods or PDE expert than I am will answer. I mostly know generalities in this area. All the Krylov methods work with evaluating Ax (and possibly A^T x) and then some operations on vectors and smaller matrices. None of the ones for nonsymmetric systems have particularly simple convergence behavior and they generally depend on finding a decent preconditioner M for which it is easy to apply M^{-1}A (or AM^{-1} for right preconditioning) and M^{-1} approximates A^{-1} in some way so that M^{-1}A has an eigenvalue distribution that improves the convergence behavior. (Although eigenvalue distribution doesn’t tell you everything for nonnormal matrices). The choice of algorithm and preconditioner can be highly dependent on the problem, which is where my own background is not sufficient. Someone here may know the right thing to do for Fokker-Planck equations, but unfortunately I am not that person.

Not having to worry about the complexity of choosing iterative methods and preconditioners for a specific problem is one argument for going with direct methods if they solve your problem with reasonable cost. Depending on the sparsity of your matrices and how good a preconditioner you can find, you might do better than the direct method, but the potential gain is much less once you have already exploited the Kronecker product structure. It might not be worth the effort, depending on your problem and your needs. Personally, I would guess that it is not likely to be worth it, but I am probably biased toward methods I’m more comfortable with and using direct methods whenever it is possible.

> [@Tamas\_Papp](#):
>
> Can you please explain or provide a link (does a package do this? [Kronecker.jl](https://juliahub.com/ui/Packages/Kronecker) seems to provide a type like this but no `\` method as far as I can tell).

This is equivalent to solving the Sylvester equation

X E\_X^T + E\_Y X = W

instead of the generalized Sylvester equation I suggested earlier. You can use `sylvc` from `MatrixEquations.jl` to solve that after forming E\_X = A\_X^{-1}C\_X and E\_Y = D\_Y^{-1} B\_Y. It would be faster than what I suggested before, but there is some potential for instability in forming W, E\_X, and E\_Y if A\_X or D\_Y is ill conditioned. The algorithm for the generalized Sylvester equation exists specifically to avoid that stability problem. If you know that A\_X and D\_Y are well conditioned, or can live with the results of moderate ill conditioning, then this would be an efficiency improvement on what I suggested. It’s a constant factor, but probably not a small constant factor. I wouldn’t be surprised if you saw a factor of 5 or 10 speedup over the generalized Sylvester equation approach.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [February 16, 2024, 5:38am UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/17 "2024-02-16T05:38:06Z")

</div>

We can’t build many preconditioners when we don’t have access to the coefficient of A.

I suggest to try `gmres` and `bicgstab` from Krylov.jl with the option `history=true`:

```julia
using Krylov

x, stats = bicgstab(A, b, history=true)

```

By default GMRES is not restarted. It’s what we want when we don’t have a good preconditioner.

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [February 16, 2024, 9:04am UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/18 "2024-02-16T09:04:17Z")

</div>

The package MatrixEquations.jl only supports Julia version\>=1.8, are there other packages also deal with generalized Sylvester equations?

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [February 16, 2024, 1:07pm UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/19 "2024-02-16T13:07:39Z")

</div>

Not that I know of. I do wonder about the reason for that lower bound. The specific code for the generalized Sylvester equation seems to have almost no dependencies. It looks like 4 functions, `gsylv`, `gsylvs!`, `gsylv2!`, `luslv!`, and `sfstruct` that don’t use anything outside of `Base` and `LinearAlgebra` and that seem likely to work for older versions of Julia. (Unless I missed something). Even if the package as a whole really does need \>=1.8 for some reason, it sure looks like if you just need to solve a generalized Sylvester equation, you could pull out those functions and use them independently of the rest of the package. It would be easy enough to check the residual to see if you got an acceptable solution. Although that’s certainly a very hacky work-around. I feel a bit bad suggesting it. And it’s even less reasonable if you want to depend on this in your own package instead of just getting the solution to an equation for your own purposes.

After writing the above, I did decide to try pulling those functions out into a separate file. The only snag in Julia 1.10.0 was that `BlasFloat` is not exported by `LinearAlgebra`.

There is also `LinearAlgebra.sylvester` if you have some well conditioned matrices so that you can convert it to an ordinary Sylvester equation. But stability will suffer if the matrices are not well conditioned.

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [February 18, 2024, 4:04am UTC](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253/20 "2024-02-18T04:04:52Z")

</div>

You had a typo here, P\_{XY}=M\_X⊗1\_Y+1\_X⊗M\_Y should be P\_{XY}=Π\_X⊗1\_Y+1\_X⊗Π\_Y

[Next page](https://discourse.julialang.org/t/solving-linear-system-with-sum-of-kronecker-products/109253.md?page=2)