# Matrix type for inversion

**URL:** https://discourse.julialang.org/t/matrix-type-for-inversion/108107
**Category:** Numerics
**Tags:** linearalgebra, sparse
**Created:** [December 27, 2023, 5:08pm UTC](https://discourse.julialang.org/t/matrix-type-for-inversion/108107 "2023-12-27T17:08:20Z")
**Posts on this page:** 3
**Page:** 1

<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: [December 27, 2023, 5:08pm UTC](https://discourse.julialang.org/t/matrix-type-for-inversion/108107/1 "2023-12-27T17:08:20Z")

</div>

Hi,

I would like to solve a linear system `A x=b` where A is band by block. The bandwith can be big (~10). An example of such matrix is for example the one in the [trapezoid](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/periodicOrbitTrapeze/) method.

I can do a LU decomposition to solve this but I am expecting that many operations will be wasted by the zero filling. Is there a way to improve this? Perhaps use a sparse matrix even if the blocks are not sparse? Use `BlockBandedMatrices` ?

Do you have any experience with this?

Thanks a lot

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [December 27, 2023, 6:21pm UTC](https://discourse.julialang.org/t/matrix-type-for-inversion/108107/2 "2023-12-27T18:21:14Z")

</div>

> [@rveltz](#):
>
> Perhaps use a sparse matrix even if the blocks are not sparse? Use `BlockBandedMatrices` ?

I would just try BandedMatrices.jl, BlockBandedMatrices.jl, and the generic SparseArrays.jl solvers, and see which ones work best for your matrix size and sparsity. (Typically, using a generic sparse solver will be slower — but only by a constant factor — than using a solver that is specialized for a banded structure, but YMMV.)

Generically, for a banded matrix (of bounded bandwidth) you can get asymptotically linear (\sim m) scaling to solve Ax=b for an m \times m matrix A, so it is definitely worth taking advantage of when m is large, compared to a dense solver which requires \sim m^3 work and \sim m^2 memory.

(A bandwidth of 10 is not particularly big, which is why I would suggest trying BandedMatrices.jl too even if the bands come in blocks.)

---

<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: [December 28, 2023, 9:22am UTC](https://discourse.julialang.org/t/matrix-type-for-inversion/108107/3 "2023-12-28T09:22:53Z")

</div>

Thank you for your answer. I will try those and see if it’s worth to write a dedicated linear solver. Basically I am wondering if it worth implementing the method of condensation of parameters (see [1,2]) to speed up the computation of periodic orbits based on a collocation method.

[1] Lust, Kurt. “Improved Numerical Floquet Multipliers.” _International Journal of Bifurcation and Chaos_ 11, no. 09 (September 2001): 2389–2410. [https://doi.org/10.1142/S0218127401003486](https://doi.org/10.1142/S0218127401003486).

[2] Doedel, Eusebius, Herbert B. Keller, and Jean Pierre Kernevez. “NUMERICAL ANALYSIS AND CONTROL OF BIFURCATION PROBLEMS (II): BIFURCATION IN INFINITE DIMENSIONS.” _International Journal of Bifurcation and Chaos_ 01, no. 04 (December 1991): 745–72. [https://doi.org/10.1142/S0218127491000555](https://doi.org/10.1142/S0218127491000555).
