Решение линейных систем уравнений с разложением SVD

Я хочу написать функцию, которая использует разложение SVD для решения системы уравнений ax=b, где a — квадратная матрица, а b — вектор значений. Функция scipy scipy.linalg.svd() должна превратить a в матрицы U W V. Для U и V я могу просто выполнить транспонирование, чтобы найти их инверсию. Но для W функция дает мне одномерный массив значений, которые мне нужно проставить по диагонали матрицы, а затем ввести единицу поверх значения.

def solveSVD(a,b):

    U,s,V=sp.svd(a,compute_uv=True)

    Ui=np.transpose(a)
    Vi=np.transpose(V)

    W=np.diag(s)

    Wi=np.empty(np.shape(W)[0],np.shape(W)[1])
    for i in range(np.shape(Wi)[0]):
        if W[i,i]!=0:
            Wi[i,i]=1/W[i,i]
    
    ai=np.matmul(Ui,np.matmul(Wi,Vi))
    x=np.matmul(ai,b)

    return(x)

Однако я получаю ошибку TypeError: тип данных не понят. Я думаю, что часть проблемы в том, что

W=np.diag(s)

не производит квадратную диагональную матрицу.

Я впервые работаю с этой библиотекой, так что извините, если я сделал что-то очень глупое, но я не могу понять, почему эта строка не работает. Спасибо всем!


python numpy scipy linear-algebra svd
person Eli Rees    schedule 11.12.2019    source источник
comment
пустая функция принимает кортеж, просто используйте W.shape, и, кроме того, ваша функция неправильно решает линейную систему   -  person Yacola    schedule 12.12.2019
comment
@Yacola Спасибо за помощь с кортежем, не могли бы вы объяснить, почему моя функция не будет правильно решать LS?   -  person Eli Rees    schedule 12.12.2019
comment
Извините, но я не уверен, как вы хотели решить эту проблему с помощью своей реализации... Является ли Ui=np.transpose(a) или более похожей на Ui=np.transpose(U), а цикл for можно просто заменить на Wi=np.diag(1/s)? Посмотрите комментарии в моем коде, чтобы понять, что делает функция.   -  person Yacola    schedule 12.12.2019

Ответы (1)


arrow_upward
4
arrow_downward

Короче говоря, использование разложения по сингулярным числам позволяет заменить исходную задачу A x = b на U diag(s) Vh x = b. Используя немного алгебры на последнем, дайте вам следующую функцию из 3 шагов, которую действительно легко читать:

import numpy as np
from scipy.linalg import svd

def solve_svd(A,b):
    # compute svd of A
    U,s,Vh = svd(A)

    # U diag(s) Vh x = b <=> diag(s) Vh x = U.T b = c
    c = np.dot(U.T,b)
    # diag(s) Vh x = c <=> Vh x = diag(1/s) c = w (trivial inversion of a diagonal matrix)
    w = np.dot(np.diag(1/s),c)
    # Vh x = w <=> x = Vh.H w (where .H stands for hermitian = conjugate transpose)
    x = np.dot(Vh.conj().T,w)
    return x

Теперь давайте проверим это с помощью

A = np.random.random((100,100))
b = np.random.random((100,1))

и сравните его с LU-разложением функции np.linalg.solve

x_svd = solve_svd(A,b)
x_lu = np.linalg.solve(A,b)

который дает

np.allclose(x_lu,x_svd)
>>> True

Пожалуйста, не стесняйтесь спрашивать дополнительные пояснения в комментариях, если это необходимо. Надеюсь это поможет.

person Yacola    schedule 11.12.2019
comment
@FatemehAsgarinejad, конструктивные комментарии всегда приветствуются, пожалуйста, не стесняйтесь редактировать этот ответ, поскольку он не предназначен быть более точным, чем solve. - person Yacola; 23.10.2020
comment
К сожалению, я не мог этого сделать. Я работаю над этим. но это дает другой ответ по сравнению со встроенной функцией для некоторых матриц. - person Fatemeh Asgarinejad; 24.10.2020