Обратная матрица numpy не работает для матрицы полного ранга - гессиан в логистической регрессии с использованием метода Ньютона

Я пытаюсь вычислить обратную матрицу полного ранга, используя numpy, но когда я тестирую точечный продукт, я обнаруживаю, что он не приводит к единичной матрице, что означает, что он не инвертировался должным образом.

Мой код:

H = calculateLogisticHessian(theta, X) #returns a 5x5 matrix
Hinv = np.linalg.inv(H)
print("H = " + str(H))
print("Hinv = " + str(Hinv))
I = np.dot(H, Hinv)
isIdentity = np.allclose(I , np.eye(5))
print("invdotinv = " + str(isIdentity) + "\n" + str(I))

и вывод:

H = [[  77.88167948   81.49914902   85.11661855   88.73408809   92.35155763]
 [  81.49914902   85.36097831   89.2228076    93.0846369    96.94646619]
 [  85.11661855   89.2228076    93.32899665   97.4351857   101.54137475]
 [  88.73408809   93.0846369    97.4351857   101.7857345   106.1362833 ]
 [  92.35155763   96.94646619  101.54137475  106.1362833   110.73119186]]
Hinv = [[  1.41918134e+02   1.00000206e+08  -1.00000632e+08  -9.99999204e+07
    1.00000205e+08]
 [  1.00000347e+08   1.00000647e+08  -4.00001421e+08   9.99994941e+07
    1.00000932e+08]
 [ -1.00000916e+08  -4.00001424e+08   8.00003700e+08   5.68436971e+02
   -3.00001928e+08]
 [ -9.99997780e+07   1.00000065e+08  -5.72321511e+02   1.00000063e+08
   -9.99997769e+07]
 [  1.00000205e+08   1.00000505e+08  -3.00001073e+08  -1.00000205e+08
    2.00000567e+08]]
invdotinv = False
[[  1.00000000e+00  -3.81469727e-06  -7.62939453e-06   3.81469727e-06
    3.81469727e-06]
 [  0.00000000e+00   1.00000191e+00  -1.52587891e-05   3.81469727e-06
    0.00000000e+00]
 [ -3.81469727e-06   1.90734863e-06   9.99992371e-01   3.81469727e-06
    3.81469727e-06]
 [  1.90734863e-06  -1.90734863e-06  -7.62939453e-06   1.00000191e+00
    3.81469727e-06]
 [  0.00000000e+00  -1.90734863e-06   0.00000000e+00   0.00000000e+00
    1.00000000e+00]]

Как видите, матрица np.dot(H, Hinv) не возвращает идентификатор и дает False при оценке np.allclose(I , np.eye(5)).

Что я делаю не так?

Позднее редактирование

это функция, которая вычисляет гессиан:

def calculateLogisticHessian(theta, X):
    '''
    calculate the hessian matrix based on given function, assuming it is some king of logistic funciton
    :param theta: the weights
    :param x: 2d array of arguments
    :return: the hessian matrix
    '''
    m, n = X.shape
    H = np.zeros((n,n))
    for i in range(0,m):
        hxi = h(theta, X[i])   #in case of logistic, will return p(y|x)
        xiDotxiT =  np.outer(X[i], np.transpose(X[i]))
        hxiTimesOneMinHxi = hxi*(1-hxi)
        currh = np.multiply(hxiTimesOneMinHxi, xiDotxiT)
        H = np.add(H, currh)
    return np.divide(H, m)

что должно быть в соответствии с формулой расчета Гессе в видео Эндрю Нга относительно метода Ньютона для логистической регрессии:

https://youtu.be/fF-6QnVB-7E?t=5m6s в 5:06

1/m * (СУММ от i=1 до m of[h(X[i]) * (1 - h(X[i]) * (X[i] * X[i]'T)])

где X — матрица данных 2x2, а h() — функция, основанная на тета (тета — вес), которая в данном случае возвращает логистическую функцию.

входы, которые я использовал:

theta = np.array([0.001, 0.002, 0.003, 0.004, 0.005])
X = np.array(range(5*7))
X = X.reshape((7,5))

H = calculateLogisticHessian(theta, X)

так есть ли ошибка в том, как я реализовал формулу Гессе, или проблема во входных данных, и в чем проблема?

Спасибо!


person Daniel Fensterheim    schedule 22.08.2017    source источник
comment
Вы уверены, что можно инвертировать эту матрицу? Детерминант приблизительно равен нулю: np.linalg.det(H) дает -5.9450505443953825e-24, что ниже числовой ненадежности, которую я бы назвал где-то около 1e-15. Также попробуйте распечатать I и прочитать документы np.allclose.   -  person Michael H.    schedule 22.08.2017
comment
Я не уверен, что матрицу можно инвертировать. Я добавил объяснение того, как я попал в эту матрицу   -  person Daniel Fensterheim    schedule 26.08.2017


Ответы (1)


Матрица Гессе часто плохо обусловлена. numpy.linalg.cond позволяет вычислить номер условия:

In [188]: np.linalg.cond(H)
Out[188]: 522295671550.72644

Поскольку число условия H велико, при вычислении его обратного значения возникают проблемы с округлением.

person B. M.    schedule 22.08.2017
comment
Я вижу, я еще не научился численному анализу. Итак, как люди обычно справляются с инверсией матрицы Гессе? Я знаю, что это требуется для таких методов, как метод Ньютона логистической регрессии. - person Daniel Fensterheim; 24.08.2017