# Solving Ax=B for large matrix dimensions efficiently in Julia

**URL:** <https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [December 9, 2020, 8:46am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504 "2020-12-09T08:46:18Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![musman](https://avatars.discourse-cdn.com/v4/letter/m/53a042/32.png) [@musman](https://discourse.julialang.org/u/musman)\
**Post date:** [December 9, 2020, 8:46am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/1 "2020-12-09T08:46:18Z")

</div>

Hi all,

I have the linear set of equation Ax=b where A can be as large as 20000 x 20000 and is invertible. What is the best computationally efficient way of solving this problem in Julia?

I have looked into linear algebra package and it has been reported that one can use A\b to get the solution? Is it computationally an efficient way to solve this problem? What about inv(A)\*b? Is there any specific package apart from linear algebra package which can handle this problem?

Looking forward for your feedbacks.

Best  
Muhammad

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [December 9, 2020, 9:02am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/2 "2020-12-09T09:02:31Z")

</div>

`x = A \ b` is the general formulation that should be efficient in most cases.

If you know that A is Hermitian and positive definite, then `x = cholesky(A) \ b` is probably better.

If you have a fast way of computing the matrix-vector product `A*x`, then have a look at [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl).

You should almost never use `x = inv(A) * b`. If you need to re-solve the problem for different right-hand sides, then it is better to compute a factorization of A than its inverse.

---

<div class="post-metadata">

**Author:** ![musman](https://avatars.discourse-cdn.com/v4/letter/m/53a042/32.png) [@musman](https://discourse.julialang.org/u/musman)\
**Post date:** [December 9, 2020, 9:09am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/3 "2020-12-09T09:09:03Z")

</div>

@Per Thanks for the reply. So, one naive question: when we use A\b, which algorithm/technique does julia use to solve it?

Yes, my equation is in a loop and I have to compute Ax=b couple of times for updated b until, the criteria is met. So, in such a case, it is better to use A\b rather than using the inv. Right?

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [December 9, 2020, 9:13am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/4 "2020-12-09T09:13:11Z")

</div>

That depends on what matrix you put in there `\` is a polyalgorithm. For square matrices julia will first perform an LU decomposition A=LU and use this decomposition to solve the problem. For rectangular `A` the result is the minimum-norm least squares solution computed by a pivoted QR factorization.

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [December 9, 2020, 9:16am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/5 "2020-12-09T09:16:14Z")

</div>

The main point beeing. If you just want to solve a single linear system of equations, using `\(A,b)` should be fine in your case. If you have multiple right hand sides `b`, then it is a good idea to precompute the factorization `Afact = factorize(A)` and use this factorization to solve your problem `\(Afact,b)`.

---

<div class="post-metadata">

**Author:** ![musman](https://avatars.discourse-cdn.com/v4/letter/m/53a042/32.png) [@musman](https://discourse.julialang.org/u/musman)\
**Post date:** [December 9, 2020, 9:22am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/6 "2020-12-09T09:22:02Z")

</div>

So, my problem is that once x is calculated (let’s say through A\b), both A and b are updated and then x is updated again based upon (A\b). The stopping criteria is when diff between x values in two consecutive iterations becomes lower than the specified threshold.

Since A and b are updated during each iteration, do you think I can still use factorize function?

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [December 9, 2020, 9:35am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/7 "2020-12-09T09:35:07Z")

</div>

Factorization is quite expensive to calculate and you would need to recalculate it in each iteration step. In this case an iterative solver as suggested by @Per would probably be faster. In general less accurate, but that is the trade off one has to make.

The only exception would be if your update to `A` analytically translates to an update of the factorization, which is usually extremely fast. An example would be if `A = M'M + aI` with identity matrix `I`, where changes of the scalar `a` have an effect on the singular values of a singular value decomposition only.

---

<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 9, 2020, 11:39am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/8 "2020-12-09T11:39:19Z")

</div>

You mean a dense matrix 20000x20000? Does it even fit in memory?

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [December 9, 2020, 11:59am UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/9 "2020-12-09T11:59:53Z")

</div>

That’s about 3 gigabytes. Not that much by modern standards, but one should take care to avoid making unnecessary copies.

---

<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 9, 2020, 1:19pm UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/10 "2020-12-09T13:19:36Z")

</div>

Where do these matrices come from? Often, large matrices in practice have some structure that you can exploit.

---

<div class="post-metadata">

**Author:** ![musman](https://avatars.discourse-cdn.com/v4/letter/m/53a042/32.png) [@musman](https://discourse.julialang.org/u/musman)\
**Post date:** [December 9, 2020, 1:23pm UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/11 "2020-12-09T13:23:36Z")

</div>

@stevengj the matrices are coming from power systems. I forgot to mention that matrix A is largely sparse. In fact, it is the jacobian matrix and has to be updated during each iteration of a while loop until a certain criteria meets.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [December 9, 2020, 2:24pm UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/12 "2020-12-09T14:24:08Z")

</div>

Is it symmetric? Positive definite? Indefinite?

---

<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 9, 2020, 2:36pm UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/13 "2020-12-09T14:36:27Z")

</div>

> [@musman](#):
>
> I forgot to mention that matrix A is largely sparse.

Then use the [sparse-matrix data structure](https://docs.julialang.org/en/v1/stdlib/SparseArrays/). `A \ b` will then use a sparse-direct solver and should be fast and memory-efficient (if `A` is sparse enough).

As mentioned above, you can save the factorization object if you are doing repeated solves with the same matrix, or use a more specialized factorization like `cholesky(A)` if your matrix has special properties like SPD.

> [@musman](#):
>
> So, my problem is that once x is calculated (let’s say through A\b), both A and b are updated and then x is updated again based upon (A\b). The stopping criteria is when diff between x values in two consecutive iterations becomes lower than the specified threshold.

It sounds like you might be solving a fixed-point equation `A(x) \ b(x) == x`? You might consider using Anderson acceleration, e.g. via the [FixedPointAcceleration package](https://github.com/s-baumann/FixedPointAcceleration.jl) or the [NLsolve package](https://github.com/JuliaNLSolvers/NLsolve.jl).

---

<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:** [January 27, 2025, 9:13pm UTC](https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimensions-efficiently-in-julia/51504/14 "2025-01-27T21:13:05Z")

</div>

2 posts were split to a new topic: [Solving AX=B with a matrix B of right-hand sides](https://discourse.julialang.org/t/solving-ax-b-with-a-matrix-b-of-right-hand-sides/125282)