<jupyter_start><jupyter_text>Why learn about LU? Why not just use `solve` (scipy) or backslash (Matlab)?It's often quite handy. For example, if you have the same matrix $A$ but many different right-hand-sides $b$, then you only need to find the LU decomposition of $A$ once (this is the expensive, $O(n^3)$ part), and then for each right-hand-side, just do backsubstitution (relatively cheap, at $O(n^2)$).For Matlab, the equivalent of `scipy.linalg.solve` is backslash `\`, aka [`mldivide`](https://www.mathworks.com/help/matlab/ref/mldivide.html).  See more at the [Matlab systems-of-equations](https://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html) documentation which also shows you how to use the new [`decomposition`](https://www.mathworks.com/help/matlab/ref/decomposition.html) function ).  You can either ask for a LU factorization explicitly and save it, or as for a vague "decomposition" and Matlab will guess a factorization (LU, Cholesky, etc) and then save it, and then you can use this with the backslash `\` seemlessly.  Neat!*To be more accurate, the true equivalent of `scipy.linalg.solve` is Matlab's [`linsolve`](https://www.mathworks.com/help/matlab/ref/linsolve.html).  It's like backslash but it doesn't guess the structure of the matrix for you, rather you have to tell it. So that can be more annoying, but it can also save computational time (if you know it's a symmetric matrix, no need for Matlab to numerically check if it is, etc.) Below is a timing example about the benefits of precomputing an LU factorization<jupyter_code>import numpy as np
from numpy.random import default_rng
rng = default_rng(123)
from timeit import timeit, time

from scipy.linalg import solve, lu_factor, lu_solve
n = int(3e3)
A = rng.standard_normal((n,n))
nTrials = 10<jupyter_output><empty_output><jupyter_text>For each new right-hand-side, just call `solve`. This is wasteful, since $A$ isn't changing<jupyter_code>%%time
for rep in range(nTrials):
  b = rng.standard_normal(n)
  x = solve(A,b)<jupyter_output>CPU times: user 12.8 s, sys: 1.42 s, total: 14.2 s
Wall time: 7.26 s<jupyter_text>This time, let's precompute the LU factorization of $A$ and re-use it. Much faster!<jupyter_code>%%time
LU, piv = lu_factor(A)
for rep in range(nTrials):
  b = rng.standard_normal(n)
  x = lu_solve((LU,piv),b)<jupyter_output>CPU times: user 1.22 s, sys: 80.9 ms, total: 1.3 s
Wall time: 682 ms
{'fraction_non_alphanumeric': 0.10839884153909805, 'fraction_numerical': 0.013239553165080678, 'mean_word_length': 4.985148514851486, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '192377', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 855, 'n_tokens_neox': 806, 'n_words': 303, 'path': 'Demos/Ch6_RepeatedSolves.ipynb', 'repo_name': 'skhadem/numerical-analysis-class', 'size': '2417'}