Крутой спуск с неоправданно большими значениями

Моя реализация наискорейшего спуска для решения Ax = b демонстрирует странное поведение: для любой достаточно большой матрицы (~ 10 x 10, пока тестировали только квадратные матрицы) возвращаемый x содержит все огромные значения (порядка 1x10^10).

def steepestDescent(A, b, numIter=100, x=None):
    """Solves Ax = b using steepest descent method"""
    warnings.filterwarnings(action="error",category=RuntimeWarning)

    # Reshape b in case it has shape (nL,)
    b = b.reshape(len(b), 1)

    exes = []
    res = []

    # Make a guess for x if none is provided
    if x==None:
        x = np.zeros((len(A[0]), 1))
        exes.append(x)

    for i in range(numIter):
        # Re-calculate r(i) using r(i) = b - Ax(i) every five iterations
        # to prevent roundoff error. Also calculates initial direction
        # of steepest descent.
        if (numIter % 5)==0:
            r = b - np.dot(A, x)
        # Otherwise use r(i+1) = r(i) - step * Ar(i)
        else:
            r = r - step * np.dot(A, r)

        res.append(r)

        # Calculate step size. Catching the runtime warning allows the function
        # to stop and return before all iterations are completed. This is
        # necessary because once the solution x has been found, r = 0, so the
        # calculation below divides by 0, turning step into "nan", which then
        # goes on to overwrite the correct answer in x with "nan"s
        try:
            step = np.dot(r.T, r) / np.dot( np.dot(r.T, A), r )
        except RuntimeWarning:
            warnings.resetwarnings()
            return x
        # Update x
        x = x + step * r
        exes.append(x)

    warnings.resetwarnings()
    return x, exes, res

(exes и res возвращаются на отладку)

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


python numpy gradient-descent numerical-methods mathematical-optimization
person Cole Zimmerman    schedule 26.07.2016    source источник
comment
Максимальный размер шага зависит от собственных значений A. Я не вижу в вашем коде вычисления собственных значений.   -  person Rodrigo de Azevedo    schedule 14.10.2016
comment
Вы не передаете step в эту функцию.   -  person kwinkunks    schedule 04.02.2018

Ответы (1)


arrow_upward
1
arrow_downward

Код кажется правильным. Например, у меня сработал следующий тест (и linalg.solve, и steepestDescent в большинстве случаев дают точный ответ):

import numpy as np

n = 100
A = np.random.random(size=(n,n)) + 10 * np.eye(n)
print(np.linalg.eig(A)[0])
b = np.random.random(size=(n,1))
x, xs, r = steepestDescent(A,b, numIter=50)
print(x - np.linalg.solve(A,b))

Проблема в математике. Этот алгоритм гарантированно сходится к правильному решению, если A - положительно определенная матрица. Добавляя единичную матрицу 10 * к случайной матрице, мы увеличиваем вероятность того, что все собственные значения положительны.

Если вы тестируете большие случайные матрицы (например, A = random.random(size=(n,n)), вы почти наверняка получите отрицательное собственное значение, и алгоритм не сойдется.

person Eolmar    schedule 28.03.2018