Понимание логики c за numpy кодом для обратного Мура-Пенроуза

Question

Я просматривал книгу под названием Практическое машинное обучение с Scikit-Learn, Keras и Tensorflow , и автор объяснял, как вычисляется псевдообратная (инверсия Мура-Пенроуза) матрицы в контексте линейной регрессии. Я цитирую здесь дословно:

Сам псевдообращение вычисляется с использованием стандартной методики факторизации матрицы, называемой декомпозицией сингулярных значений (SVD), которая может разложить матрицу обучающего набора X на умножение матриц трех матриц U Σ V T (см. numpy .linalg.svd ()). Псевдообращение вычисляется как X + = V * Σ + * U T. Чтобы вычислить матрицу Σ +, алгоритм берет Σ и устанавливает в ноль все значения, меньшие порогового значения, а затем заменяет все ненулевые значения на их обратные и, наконец, транспонирует Полученная матрица. Этот подход более эффективен, чем вычисление уравнения Нормы.

Я получил представление о том, как псевдообратная связь и SVD связаны с этой публикацией. Но я не могу найти asp обоснование для установки всех значений меньше порогового значения в ноль. Обратная диагональная матрица получается путем взятия обратных значений диагональных элементов. Тогда маленькие значения будут преобразованы в большие значения в обратной матрице, верно? Тогда почему мы удаляем большие значения?

Я пошел и посмотрел код numpy, и он выглядит следующим образом, просто для справки:

@array_function_dispatch(_pinv_dispatcher)
def pinv(a, rcond=1e-15, hermitian=False):
a, wrap = _makearray(a)
    rcond = asarray(rcond)
    if _is_empty_2d(a):
        m, n = a.shape[-2:]
        res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
        return wrap(res)
    a = a.conjugate()
    u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)

    # discard small singular values
    cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
    large = s > cutoff
    s = divide(1, s, where=large, out=s)
    s[~large] = 0

    res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
    return wrap(res)

Answer 1

Это почти наверняка поправка на числовую ошибку. Чтобы понять, почему это может быть необходимо, посмотрите, что происходит, когда вы берете svd матрицы 2x2 ранга один. Мы можем создать матрицу ранга один, взяв внешнее произведение вектора следующим образом:

>>> a = numpy.arange(2) + 1
>>> A = a[:, None] * a[None, :]
>>> A
array([[1, 2],
       [2, 4]])

Несмотря на то, что это матрица 2x2, она имеет только один линейно независимый столбец, поэтому ее ранг один вместо двух. Поэтому следует ожидать, что когда мы передадим его в svd, одно из сингулярных значений будет равно нулю. Но посмотрите, что происходит:

>>> U, s, V = numpy.linalg.svd(A)
>>> s
array([5.00000000e+00, 1.98602732e-16])

На самом деле мы получаем единственное значение, которое не совсем ноль. Этот результат неизбежен во многих случаях, учитывая, что мы работаем с числами с плавающей запятой конечной точности. Таким образом, хотя проблема, которую вы определили, является реальной, мы не сможем на практике определить разницу между матрицей, которая действительно имеет очень маленькое сингулярное значение, и матрицей, которая должна иметь нулевое сингулярное значение, но не имеет. Установка небольших значений на ноль - самый безопасный практический способ решения этой проблемы.