

In [1]:

    
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
import time

%matplotlib inline



In [2]:

    
np.random.seed(17)

Triagonal matrices



In [3]:

    
def generate_random_triagonal(n):
    A = np.zeros((n, n))
    for i in range(n - 1):
        A[i + 1][i] = rand()
        A[i][i + 1] = rand()
        A[i][i] = rand()
    A[n - 1][n - 1] = rand()
    return A, rand(n)



In [4]:

    
def norm(A, x, f):
    return np.linalg.norm(A.dot(x) - f)



In [5]:

    
def run(method, A, f, verbose=False):
    if not verbose:
        print("-" * 100)
        print(method.__name__.upper())
    x = method(A, f)
    score = norm(A, x, f)
    if not verbose:
        print("x =", x)
        print("score =", score)
    return score



In [6]:

    
def linalg(A, f):
    return np.linalg.solve(A, f)



In [7]:

    
def right(A, f):    
    n = A.shape[0]
    a = np.zeros(n - 1)
    b = np.zeros(n - 1)
    a[0] = -A[0][1] / A[0][0]
    b[0] = f[0] / A[0][0]
    for i in range(1, n - 1):
        a[i] = -A[i][i + 1] / (A[i][i - 1] * a[i - 1] + A[i][i])
        b[i] = (f[i] - A[i][i - 1] * b[i - 1]) / (A[i][i - 1] * a[i - 1] + A[i][i])
    x = np.zeros(n)
    x[n - 1] = (f[n - 1] - A[n - 1][n - 2] * b[n - 2]) / (A[n - 1][n - 2] * a[n - 2] + A[n - 1][n - 1])
    for i in range(n - 2, -1, -1):
        x[i] = x[i + 1] * a[i] + b[i]
    return x



In [8]:

    
def left(A, f):    
    n = A.shape[0]
    a = np.zeros(n)
    b = np.zeros(n)
    a[n - 1] = -A[n - 1][n - 2] / A[n - 1][n - 1]
    b[n - 1] = f[n - 1] / A[n - 1][n - 1]
    for i in range(n - 2, 0, -1):
        a[i] = -A[i][i - 1] / (A[i][i + 1] * a[i + 1] + A[i][i])
        b[i] = (f[i] - A[i][i + 1] * b[i + 1]) / (A[i][i + 1] * a[i + 1] + A[i][i])
    x = np.zeros(n)
    x[0] = (f[0] - A[0][1] * b[1]) / (A[0][1] * a[1] + A[0][0])
    for i in range(1, n):
        x[i] = x[i - 1] * a[i] + b[i]
    return x

2x2 corner case



In [9]:

    
A = np.array(
    [[1, 2],
     [3, 4]]
)
f = np.array([5, 11])

_ = run(right, A, f)
_ = run(left, A, f)









    



----------------------------------------------------------------------------------------------------
RIGHT
x = [ 1.  2.]
score = 0.0
----------------------------------------------------------------------------------------------------
LEFT
x = [ 1.  2.]
score = 0.0

Another example



In [10]:

    
A, f = generate_random_triagonal(6)
print(A)

_ = run(right, A, f)
_ = run(left, A, f)









    



[[ 0.19152079  0.53058676  0.          0.          0.          0.        ]
 [ 0.294665    0.65633352  0.78698546  0.          0.          0.        ]
 [ 0.          0.06790036  0.03906292  0.57560289  0.          0.        ]
 [ 0.          0.          0.6375209   0.06004468  0.94568319  0.        ]
 [ 0.          0.          0.          0.3578136   0.05119367  0.87729053]
 [ 0.          0.          0.          0.          0.8640421   0.65241862]]
----------------------------------------------------------------------------------------------------
RIGHT
x = [ 10.25298105  -2.66103067  -0.86044513   1.21233764   0.80232523
  -0.20191618]
score = 1.84109660315e-16
----------------------------------------------------------------------------------------------------
LEFT
x = [ 10.25298105  -2.66103067  -0.86044513   1.21233764   0.80232523
  -0.20191618]
score = 1.46868701149e-16

Huge matrices



In [11]:

    
def timeit(method, A, f):
    start = time.process_time()
    method(A, f)
    return time.process_time() - start



In [12]:

    
A, f = generate_random_triagonal(5000)

%timeit run(right, A, f, verbose=True)
%timeit run(left, A, f, verbose=True)
%timeit run(linalg, A, f, verbose=True)









    



10 loops, best of 3: 40.5 ms per loop
10 loops, best of 3: 49.6 ms per loop
1 loop, best of 3: 1.31 s per loop



In [13]:

    
def plot_speed(method, sizes):
    scores = np.zeros_like(sizes, dtype=np.float)
    for i in range(len(sizes)):
        A, f = generate_random_triagonal(sizes[i])
        scores[i] = timeit(method, A, f)
    plt.plot(sizes, scores, label=method.__name__)



In [14]:

    
sizes = np.linspace(2, 5000, num=50, dtype=np.int)

plt.figure(figsize=(15, 10))
plot_speed(right, sizes)
plot_speed(left, sizes)
plot_speed(linalg, sizes)
plt.title("Time comparison of Right/Left/Linalg methods for random triagonal matrices").set_fontsize("xx-large")
plt.xlabel("size").set_fontsize("xx-large")
plt.ylabel("time").set_fontsize("xx-large")
# plt.xscale("log")
legend = plt.legend(loc="upper right")
for label in legend.get_texts():
    label.set_fontsize("xx-large")
plt.show()