Least Squares Minimal Norm Solution && SVD

 Начнем с хорошего примера. Чтобы иметь контроль над матрицей, мы строим ее с помощью ее разложения по сингулярным числам  (SVD) A = USV' с ортогональными матрицами U и V и диагональной матрицей S. Напомним, что S имеет только неотрицательные элементы и — для правильных матриц — даже строго положительные значения. Поэтому мы начинаем с установки 

(S) = (1, 2, 3, 4, 5).

Помимо этого приведен пример работы np.linalg.lstsq с другими данными

(.env) boris@boris-All-Series:~/MATPLOTLIBSR/SVD$ cat scipySVD1.py

from pprint import pprint

import matplotlib.pyplot as plt

import numpy as np

from numpy.linalg import norm

from scipy.stats import ortho_group

from scipy import linalg

from scipy.sparse import linalg as spla

def generate_U_S_Vt(n=10, p=5, random_state=532):

    """

    Сгенерируйте SVD, чтобы построить правильную матрицу A.

    A имеет n строк, p столбцов.

    Возвращает

    U: ортогональная матрица

    S: диагональная матрица

    Vt: ортогональная матрица

    """     

    r = min(n, p)

    S = np.diag(1.0 * np.arange(1, 1 + r))

    if n > p:

        # add rows with value 0

        S = np.concatenate((S, np.zeros((n - p, p))), axis=0)

    elif p > n:

        # add columns with value 0

        S = np.concatenate((S, np.zeros((n, p - n))), axis=1)

    U = ortho_group.rvs(n, random_state=random_state)

    Vt = ortho_group.rvs(p, random_state=random_state + 1)

    return U, S, Vt

def solve_least_squares(A, b):

    """

    Решите метод наименьших квадратов несколькими методами.

    Возвращает

    x : словарь с решателем и решением

    """

    x = {}

    x["gelsd"] = linalg.lstsq(A, b, lapack_driver="gelsd")[0]

    x["gelsy"] = linalg.lstsq(A, b, lapack_driver="gelsy")[0]

    x["lsqr"] = spla.lsqr(A, b)[0]

    x["lsmr"] = spla.lsmr(A, b)[0]

    x["normal_eq"] = linalg.solve(A.T @ A, A.T @ b, assume_a="sym")

    return x

def print_dict(d):

    np.set_string_function(np.array2string)

    pprint(d)

    np.set_string_function(None)

np.set_printoptions(precision=5)

n = 10

p = 5

U, S, Vt = generate_U_S_Vt(n=n, p=p)

A = U @ S @ Vt

x_true = np.round(6 * Vt.T[:p, 0])  # interesting choice

rng = np.random.default_rng(157)

noise = rng.standard_normal(n)

b = A @ x_true + noise

S_inv = np.copy(S.T)

S_inv[S_inv>0] = 1/S_inv[S_inv>0]

x_exact = Vt.T @ S_inv @ U.T @ b

print(f"x_exact = {x_exact}")

print("\n")

print_dict(solve_least_squares(A, b))

(.env) boris@boris-All-Series:~/MATPLOTLIBSR/SVD$ python3 scipySVD1.py

x_exact = [ 0.78087 -4.74942 -0.99938 -2.38327 -3.7431 ]

{'gelsd': [ 0.78087 -4.74942 -0.99938 -2.38327 -3.7431 ],

 'gelsy': [ 0.78087 -4.74942 -0.99938 -2.38327 -3.7431 ],

 'lsmr': [ 0.78087 -4.74942 -0.99938 -2.38327 -3.7431 ],

 'lsqr': [ 0.78087 -4.74942 -0.99938 -2.38327 -3.7431 ],

 'normal_eq': [ 0.78087 -4.74942 -0.99938 -2.38327 -3.7431 ]}

Программа для демонстрации работы np.linalg.lstsq

(.env) boris@boris-All-Series:~/MATPLOTLIBSR/SVD$ cat  npLinalgLstsq.py

import numpy as np

import matplotlib.pyplot as plt

import warnings

warnings.filterwarnings("ignore")

A = np.array([[1, 2, 1],

              [1,1,2],

              [2,1,1],

              [1,1,1]])

b = np.array([4,3,5,4])

x, residuals, rank, s = np.linalg.lstsq(A,b)

print(x)

(.env) boris@boris-All-Series:~/MATPLOTLIBSR/SVD$ python3  npLinalgLstsq.py

[2.05263158 1.05263158 0.05263158]

(.env) [boris@sever35fedora SVD]$ cat npLinalgLstsq.py

import numpy as np

import matplotlib.pyplot as plt

A = np.array([[1, 2, 1],

              [1,1,2],

              [2,1,1],

              [1,1,1]])

b = np.array([4,3,5,4])

x, residuals, rank, s = np.linalg.lstsq(A,b,rcond=-1)

print(x)

(.env) [boris@sever35fedora SVD]$ python npLinalgLstsq.py

[2.05263158 1.05263158 0.05263158]

(.env) boris@boris-All-Series:~/MATPLOTLIBSR/SVD$ cat npLinalgLstsq1.py

import numpy as np

import matplotlib.pyplot as plt

A = np.array([[1, 2, 1],

              [1,1,2],

              [2,1,1],

              [1,1,1]])

# Eсли b имеет более одного измерения, 

# lstsq решит систему, соответствующую каждому столбцу b

 

b = np.array([[4,3,5,4],[1,2,3,4]]).T

x,residuals,rank,s = np.linalg.lstsq(A,b,rcond=-1)

print(x)

(.env) boris@boris-All-Series:~/MATPLOTLIBSR/SVD$ python3  npLinalgLstsq1.py

[[ 2.05263158  1.63157895]

 [ 1.05263158 -0.36842105]

 [ 0.05263158  0.63157895]]