# 9.3. Linear Systems and Regression#

using LinearAlgebra


## 9.3.1. Linear Systems#

One of the most common uses of matrices is for solving linear systems of equations. Julia uses the backslash operator \ for this:

A = [1 2; 3 4]
b = [5,1]
x = A \ b         # Solve Ax = b for x

2-element Vector{Float64}:
-9.0
7.0


One way to view the syntax A\b is that it multiplies by A-inverse from the left, but using much more efficient and accurate algorithms.

To check that the answer is correct, we can use the $$\approx$$ operator (since floating point numbers should never be compared exactly):

A*x ≈ b          # Confirm solution is correct

true


For systems with many right-hand side vectors b, the \ operator also works with matrices:

B = [5 7; 1 -3]
X = A \ B        # Solve for two RHS vectors

2×2 Matrix{Float64}:
-9.0  -17.0
7.0   12.0

 A*X ≈ B         # Confirm solution is correct

true


The algorithm used by the \ operator is typically Gaussian elimination, but the details are quite complex depending on the type of matrices involved. Due to the high cost of general Gaussian elimination, it can make a big difference if you use a specialized matrix type:

n = 2000
T = SymTridiagonal(2ones(n), -ones(n))     # n-by-n symmetric tridiagonal

for rep = 1:3 @time T \ randn(n) end       # Very fast since T is a SymTridiagonal
Tfull = Matrix(T)                          # Convert T to a full 2D array
for rep = 1:3 @time Tfull \ randn(n) end   # Now \ is magnitudes slower

  0.134713 seconds (389.64 k allocations: 23.538 MiB, 99.96% compilation time)
0.000025 seconds (4 allocations: 63.000 KiB)
0.000017 seconds (4 allocations: 63.000 KiB)
0.216786 seconds (295.09 k allocations: 50.449 MiB, 3.38% gc time, 57.38% compilation time)
0.104367 seconds (5 allocations: 30.564 MiB, 12.76% gc time)
0.090345 seconds (5 allocations: 30.564 MiB)


The matrix A in A\b can also be rectangular, in which case a minimum-norm least squares solution is computed.

## 9.3.2. Linear regression#

Suppose you want to approximate a set of $$n$$ points $$(x_i,y_i)$$, $$i=1,\ldots,n$$, by a straight line. The least squares approximation $$y=a + bx$$ is given by the least-squares solution of the following over-determined system:

$\begin{split} \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix}= \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} \end{split}$
x = 0:0.1:10
n = length(x)
y = 3x .- 2 + randn(n)     # Example data: straight line with noise

A = [ones(n) x]            # LHS
ab = A \ y                 # Least-squares solution

using PyPlot
xplot = 0:10;
yplot = @. ab[1] + ab[2] * xplot
plot(x,y,".")
plot(xplot, yplot, "r");