linsys_solve_demo2.ipynb

Kernel: Python 3 (ipykernel)

#  Solving  linear systems


import numpy as np
from scipy import linalg
from itertools import combinations
import scipy


def naive_solve(A, b):
    return np.linalg.inv(A.T @ A) @ A.T @ b


def qr_solve(A, b):
    m, n = np.shape(A)
    if m > n:
        return qr_solve_over(A, b)
    else:
        return qr_solve_under(A, b)


def qr_solve_over(A, b):
    Q, R = np.linalg.qr(A)
    Qb = np.dot(Q.T, b)
    return scipy.linalg.solve_triangular(R, Qb)


def qr_solve_under(A, b):
    Q, R = np.linalg.qr(A.T)
    c = scipy.linalg.solve_triangular(R, b)
    m, n = A.shape
    x = np.zeros(n)
    xx = np.dot(Q, c)
    xx = xx[:, 0]
    K = xx.shape[0]
    x[:K] = xx  # other components are zero
    return x


def run_demo(m, n):
    print("Solving linear system with {} constraints and {} unknowns".format(m, n))
    np.random.seed(0)
    A = np.random.rand(m, n)
    b = np.random.rand(m, 1)

    methods = list()
    solns = list()

    methods.append("naive")
    solns.append(naive_solve(A, b))

    methods.append("pinv")
    solns.append(np.dot(np.linalg.pinv(A), b))

    methods.append("lstsq")
    solns.append(np.linalg.lstsq(A, b, rcond=None)[0])

    methods.append("qr")
    solns.append(qr_solve(A, b))

    for (method, soln) in zip(methods, solns):
        residual = b - np.dot(A, soln)
        print(
            "method {}, norm {:0.5f}, residual {:0.5f}".format(method, np.linalg.norm(soln), np.linalg.norm(residual))
        )
        print(soln.T)

    # https://stackoverflow.com/questions/33559946/numpy-vs-mldivide-matlab-operator
    if m < n:  # underdetermined
        rank = np.linalg.matrix_rank(A)
        assert m == rank
        for nz in combinations(range(n), rank):  # the variables not set to zero
            soln = np.zeros((n, 1))
            soln[nz, :] = np.asarray(np.linalg.solve(A[:, nz], b))
            residual = b - np.dot(A, soln)
            print("sparse qr, norm {:0.5f}, residual {:0.5f}".format(np.linalg.norm(soln), np.linalg.norm(residual)))
            print(soln.T)


m = 5
n = 3  # Overdetermined
run_demo(m, n)

m = 3
n = 5  # Underdetermined
run_demo(m, n)

Solving linear system with 5 constraints and 3 unknowns
method naive, norm 2.20278, residual 0.33863
[[-1.37054293  1.7166056  -0.1646944 ]]
method pinv, norm 2.20278, residual 0.33863
[[-1.37054293  1.7166056  -0.1646944 ]]
method lstsq, norm 2.20278, residual 0.33863
[[-1.37054293  1.7166056  -0.1646944 ]]
method qr, norm 2.20278, residual 0.33863
[[-1.37054293  1.7166056  -0.1646944 ]]
Solving linear system with 3 constraints and 5 unknowns
method naive, norm 2.35511, residual 0.39906
[[ 0.3265749   1.13083441 -0.38484624  0.63303706 -1.90059589]]
method pinv, norm 1.85527, residual 0.00000
[[ 0.95039618  0.48809878 -0.86619904  0.42900372 -1.16884755]]
method lstsq, norm 1.85527, residual 0.00000
[[ 0.95039618  0.48809878 -0.86619904  0.42900372 -1.16884755]]
method qr, norm 2.41014, residual 1.66621
[ 1.09891785  1.55804808 -1.18142201 -0.16370947 -0.86664275]
sparse qr, norm 2.95210, residual 0.00000
[[ 2.47618906 -0.48751645 -1.53156305  0.          0.        ]]
sparse qr, norm 7.27349, residual 0.00000
[[ 6.18106604 -2.26435156  0.         -3.09365953  0.        ]]
sparse qr, norm 2.29785, residual 0.00000
[[ 0.71032114  0.78487888  0.          0.         -2.03949434]]
sparse qr, norm 2.58081, residual 0.00000
[[ 1.45966904  0.         -1.95178335  0.84881816  0.        ]]
sparse qr, norm 2.17755, residual 0.00000
[[ 1.79959929  0.         -0.94474686  0.         -0.78142934]]
sparse qr, norm 2.72323, residual 0.00000
[[ 2.11850337  0.          0.         -0.79631503 -1.51452383]]
sparse qr, norm 3.36071, residual 0.00000
[[ 0.          0.70004786 -2.55519749  2.06767619  0.        ]]
sparse qr, norm 3.19998, residual 0.00000
[[ 0.          1.29670046  0.6160719   0.         -2.85988182]]
sparse qr, norm 2.62020, residual 0.00000
[[ 0.          1.18079075  0.          0.40168054 -2.30430214]]
sparse qr, norm 8.41669, residual 0.00000
[[ 0.          0.         -6.27602323  4.49366459  3.35547712]]