# Using Newton Method in a Neural Network/

I programmed a Neural Network to do binary classification in python, and during the backpropagation step I used Newton-Raphson's method for optimization. Any kind of feedback would be appreciated, but mostly I'd like to know if all of the gradients and hessians were computed correctly - for the examples I ran, it seems to be working correctly, leading me to believe everything is ok. Another major help I needed is on how to make the code more efficient - the computational cost for creating the hessian matrices is really high as it is, I wish there was a faster way of doing it.

An observation on the hessian: to guarantee that the update direction for the weights and biases is a descent direction, I used aprox_pos_def() function to approximate the hessian with a positive definite matrix; you can see the theory behind this transformation here http://www.optimization-online.org/DB_FILE/2003/12/800.pdf.

For a better understanding of the Newton-Raphson method, these wikipedia articles are good enough:

Here is the code:

#!/usr/bin/env python3
# Author: Erik Davino Vincent

# Almost pure Newton-Raphson method implementation for a neural network

import numpy as np
import matplotlib.pyplot as plt

# Fits for X and Y
def model(X, Y, layerSizes, learningRate, max_iter = 100, plotN = 100):

n_x, m = X.shape
n_y = Y.shape
if m != Y.shape:
raise ValueError("Invalid vector sizes for X and Y -> X size = " + str(X.shape) + " while Y size = " + str(Y.shape) + ".")

weights = initializeWeights(n_x, n_y, layerSizes)
numLayers = len(layerSizes)

return newtonMethod(X, Y, weights, learningRate, numLayers, max_iter, plotN)

# Sigmoid activation function
def sigmoid(t):
return 1/(1+np.exp(-t))

# Initializes weights for each layer
def initializeWeights(n_x, n_y, layerSizes, scaler = 0.1, seed = None):

np.random.seed(seed)

weights = {}
n_hPrev = n_x
for i in range(len(layerSizes)):
n_h = layerSizes[i]
weights['W'+str(i+1)] = np.random.randn(n_h, n_hPrev)*scaler
weights['b'+str(i+1)] = np.zeros((n_h,1))
n_hPrev = n_h

weights['W'+str(len(layerSizes)+1)] = np.random.random((n_y, n_hPrev))*scaler
weights['b'+str(len(layerSizes)+1)] = np.zeros((n_y,1))

np.random.seed(None)

return weights

# Creates a positive definite aproximation to a not positive definite matrix
def aprox_pos_def(A):

u, V = np.linalg.eig(A)
U = np.absolute(np.diag(u))
B = np.dot(V, np.dot(U, V.T))

return B

# Performs foward propagation
def forwardPropagation(weights, X, numLayers):

# Cache for A
Avals = {}
AiPrev = np.copy(X)
for i in range(numLayers+1):

Wi = weights['W'+str(i+1)]
bi = weights['b'+str(i+1)]

Zi = np.dot(Wi, AiPrev) + bi
Ai = sigmoid(Zi)

Avals['A'+str(i+1)] = Ai
AiPrev = np.copy(Ai)

return Avals

# Newton method for training a neural network
def newtonMethod(X, Y, weights, learningRate, numLayers, max_iter, plotN):

# Cache for ploting cost
cost = []
iteration = []

# Cache for minimum cost and best weights
bestWeights = weights.copy()
minCost = np.inf

# Init break_code = 0
break_code = 0
for it in range(max_iter):

# Forward propagation
Avals = forwardPropagation(weights, X, numLayers)

# Evaluates cost fucntion
AL = np.copy(Avals['A'+str(numLayers+1)])
lossFunc = -(Y*np.log(AL) + (1-Y)*np.log(1-AL))
costFunc = np.mean(lossFunc)

if costFunc < minCost:
bestWeights = weights.copy()
minCost = costFunc

# Caches cost function every plotN iterations
if it%plotN == 0:
cost.append(costFunc)
iteration.append(it)

# "Backward" propagation (Newton-Raphson's Method) loop
for i in range(numLayers+1, 0, -1):

# Gets current layer weights
Wi = np.copy(weights['W'+str(i)])
bi = np.copy(weights['b'+str(i)])
Ai = np.copy(Avals['A'+str(i)])

# If on the first layer, APrev = X; else APrev = Ai-1
if i == 1:
APrev = np.copy(X)
else:
APrev = np.copy(Avals['A'+str(i-1)])

# If on the last layer, dZi = Ai - Y; else dZi = (Wi+1 . dZi+1) * (Ai*(1-Ai))
if i == numLayers+1:
dZi = (Ai - Y)  # /(Ai * (1 - Ai)) ???
else:
dZi = np.dot(Wnxt.T, dZnxt) * Ai * (1 - Ai)

# Calculates gradient vector (actually a matrix) of i-th layer
m = APrev.shape
gradbi = np.sum(dZi, axis = 1, keepdims = 1)/m

# Performs newton method on each node of i-th layer
try:
for j in range(len(Ai)):

# Creates hessian matrix for node j in layer i
hessMatxi = np.dot(np.array([Ai[j]]), (1-np.array([Ai[j]])).T) * np.dot(APrev, APrev.T)/m
hessbipar = np.dot(np.array([Ai[j]]), (1-np.array([Ai[j]])).T) * np.dot(APrev, np.ones((APrev.shape,1)))/m
hessbi = np.dot(np.array([Ai[j]]), (1-np.array([Ai[j]])).T)/m
hessMatxi = np.concatenate((hessMatxi, hessbipar), axis = 1)
hessbipar = np.concatenate((hessbipar, hessbi), axis = 0)
hessMatxi = np.concatenate((hessMatxi, hessbipar.T), axis = 0)
hessMatxi = aprox_pos_def(hessMatxi)

# Creates descent direction for layer i
if j == 0:
else:
deltai = np.append(deltai, np.linalg.solve(hessMatxi, np.array([gradVecti[:,j]]).T).T, axis = 0)
except:
print("Singular matrix found when calculating descent direction; terminating computation.")
break_code = 1
break

# Descent step for weights and biases
dWi = deltai[:,:-1]
dbi = np.array([deltai[:,-1]]).T

# Cache dZi, Wi, bi
dZnxt = np.copy(dZi)
Wnxt = np.copy(Wi)

Wi = Wi - learningRate*dWi
bi = bi - learningRate*dbi
weights['W'+str(i)] = Wi
weights['b'+str(i)] = bi

# Plot cost every plotN iterations
if it % plotN == 0:
plt.clf()
plt.plot(iteration, cost, color = 'b')
plt.xlabel("Iteration")
plt.ylabel("Cost Function")
plt.title("Newton-Raphson Descent for Cost Function")
plt.pause(0.001)

# End of backprop loop ==================================================

# Early breaking condition met
if break_code:
AL = Avals['A'+str(numLayers+1)]
lossFunc = -(Y*np.log(AL) + (1-Y)*np.log(1-AL))
costFunc = np.mean(lossFunc)
cost.append(costFunc)
iteration.append(it)
break

# Plots cost function over time
plt.clf()
plt.plot(iteration, cost, color = 'b')
plt.xlabel("Iteration")
plt.ylabel("Cost Function")
plt.title("Newton-Raphson Descent for Cost Function")
plt.show(block = 0)

return bestWeights

# Predicts if X vector tag is 1 or 0
def predict(weights, X, numLayers):

A = forwardPropagation(weights, X, numLayers)
A = A['A'+str(numLayers+1)]

return (A > 0.5)

def example():

from sklearn.preprocessing import StandardScaler as StdScaler
from sklearn.model_selection import train_test_split

# Won't work without scalling data
scaler = StdScaler()

X = data
scaler.fit(X)
X = scaler.transform(X)
y = data

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.30)

X_train = X_train.T
X_test = X_test.T
y_train = np.array([y_train])
y_test = np.array([y_test])

layers = [5, 5, 5]
weights = model(X_train, y_train, layers, 0.1, max_iter = 500, plotN = 10)

pred = predict(weights, X_train, len(layers))
percnt = 0
for i in range(pred.shape):
if pred[0,i] == y_train[0,i]:
percnt += 1
percnt /= pred.shape
print()
print("Accuracy of", percnt*100 , "% on training set")

pred = predict(weights, X_test, len(layers))
percnt = 0
for i in range(pred.shape):
if pred[0,i] == y_test[0,i]:
percnt += 1
percnt /= pred.shape
print()
print("Accuracy of", percnt*100 , "% on test set")

example()


## Consistent unpacking

Rather than sometimes using indexes, can you do

n_x, m_x = X.shape
n_y, m_y = Y.shape


?

## Interpolation

"Invalid vector sizes for X and Y -> X size = " + str(X.shape) + " while Y size = " + str(Y.shape) + "."


can be

f'Invalid vector sizes for X and Y -> X size = {X.shape} while Y size = {Y.shape}.'


This:

"Accuracy of", percnt*100 , "% on training set"


should similarly use an f-string, and does not need to multiply by 100 if you use the built-in percent field type.

## snake_case

By convention,

initializeWeights


should be

initialize_weights


and similar for your other functions and local variables. In particular Avals should be lower-case; otherwise it looks like class.

## Expressions

Expressions such as

dZi = (Ai - Y)
return (A > 0.5)


do not need parens.

## Bare except:

This is dangerous and prevents user break (Ctrl+C) from working. Instead, except Exception, or ideally something more specific if you know what you expect to see.

## In-place subtraction

Wi = Wi - learningRate*dWi
bi = bi - learningRate*dbi


can use -=. This will improve brevity and may marginally improve performance.

## Equivalent equations

It's minor, but

-(Y*np.log(AL) + (1-Y)*np.log(1-AL))


is equivalent to

(Y - 1)*np.log(1 - AL) - Y*np.log(AL)