Here are some of the steps that need to happen before a release:

• Add tests for everything
• Move factorizations to CySparse
• Move iterative methods to PyKrylov
• Move NLPy to the optimizers organization
• Integrate the simplified slack model
• Use try blocks to check for optional dependencies at run time
• precon/pycfs will be replaced by lldl
• ensure all code is PEP8 and PEP257 compliant

@syarra @counterclocker

nlpy/model tests don't run anymore. It seems nose is trying to import _amplpy locally but isn't searching PYTHONPATH. The module works from the command line.

Perhaps we shouldn't use numpy.testing, which only allows to do import nlpy; nlpy.test().

Add continuous integration with Travis and (perhaps) Appveyor. Have nosetests report code coverage.

import numpy as np
from keras.datasets import mnist
from sklearn.model_selection import train_test_split

Load MNIST dataset

(X, y), (X_test, y_test) = mnist.load_data()

Subset data to use only class 0 and class 1

indices = np.logical_or(y == 0, y == 1)
X = X[indices]
y = y[indices]

Reshape images to 1D vectors

X = X.reshape(X.shape[0], -1)

Standardize dataset

mean = np.mean(X)
std = np.std(X)
X_std = (X - mean) / std

Split data into training and validation sets

X_train, X_val, y_train, y_val = train_test_split(X_std, y, test_size=0.2, random_state=42)

Define the sigmoid function

def sigmoid(z):
return 1 / (1 + np.exp(-z))
learning_rate = 0.01
num_iterations = 1000

Use L1 regularization with gradient descent optimizer:

lambdas = [0.001, 0.01]

for l in lambdas:
# Initialize the parameters
w = np.zeros(X_train.shape[1])
b = 0

for i in range(num_iterations):
    # Forward pass
    z = np.dot(X_train, w) + b
    y_pred = sigmoid(z)

    # Compute the cost
    cost = -np.mean(y_train * np.log(y_pred) + (1 - y_train) * np.log(1 - y_pred)) + l * np.sum(np.abs(w))

    # Backward pass
    dz = y_pred - y_train
    dw = np.dot(X_train.T, dz) / X_train.shape[0] + l * np.sign(w)
    db = np.mean(dz)

    # Update the parameters
    w = w - learning_rate * dw
    b = b - learning_rate * db

# Evaluate the model on the validation data
z_val = np.dot(X_val, w) + b
y_val_pred = sigmoid(z_val)
y_val_pred[y_val_pred >= 0.5] = 1
y_val_pred[y_val_pred < 0.5] = 0
accuracy = np.mean(y_val_pred == y_val)
print("Lambda: {}, Validation accuracy: {}".format(l, accuracy))

Use mini-batch gradient descent optimizer:

batch_sizes = [128, 64]

Define the mini-batch generator function

def minibatch_generator(X, y, batch_size):
num_samples = X.shape[0]
indices = np.arange(num_samples)
np.random.shuffle(indices)
for start_idx in range(0, num_samples - batch_size + 1, batch_size):
excerpt = indices[start_idx:start_idx + batch_size]
yield X[excerpt], y[excerpt]

Train the logistic regression model using mini-batch gradient descent

for batch_size in batch_sizes:
print(f"Batch size: {batch_size}")
w = np.zeros(X_train.shape[1])
b = 0
for i in range(num_iterations):
# Mini-batch generator
batch_generator = minibatch_generator(X_train, y_train, batch_size)

    for batch_X, batch_y in batch_generator:
        # Forward pass
        z = np.dot(batch_X, w) + b
        y_pred = sigmoid(z)

        # Compute the cost
        cost = -np.mean(batch_y * np.log(y_pred) + (1 - batch_y) * np.log(1 - y_pred))

        # Backward pass
        dz = y_pred - batch_y
        dw = np.dot(batch_X.T, dz) / batch_size
        db = np.mean(dz)

        # Update the parameters
        w = w - learning_rate * dw
        b = b - learning_rate * db

# Evaluate the model on the validation data
z_val = np.dot(X_val, w) + b
y_val_pred = sigmoid(z_val)
y_val_pred[y_val_pred >= 0.5] = 1
y_val_pred[y_val_pred < 0.5] = 0
accuracy = np.mean(y_val_pred == y_val)
print("Validation accuracy:", accuracy)

#RMS Prob optimizer
eps = 1e-8
beta = 0.9
s_w = np.zeros(X_train.shape[1])
s_b = 0

for i in range(num_iterations):
# Forward pass
z = np.dot(X_train, w) + b
y_pred = sigmoid(z)

# Compute the cost
cost = -np.mean(y_train * np.log(y_pred) + (1 - y_train) * np.log(1 - y_pred))

# Backward pass
dz = y_pred - y_train
dw = np.dot(X_train.T, dz) / X_train.shape[0]
db = np.mean(dz)

# Update the RMSprop parameters
s_w = beta * s_w + (1 - beta) * np.square(dw)
s_b = beta * s_b + (1 - beta) * np.square(db)

# Update the parameters using RMSprop optimizer
w = w - learning_rate * dw / np.sqrt(s_w + eps)
b = b - learning_rate * db / np.sqrt(s_b + eps)

Evaluate the model on the validation data

z_val = np.dot(X_val, w) + b
y_val_pred = sigmoid(z_val)
y_val_pred[y_val_pred >= 0.5] = 1
y_val_pred[y_val_pred < 0.5] = 0
accuracy = np.mean(y_val_pred == y_val)
print("Validation accuracy:", accuracy)

Adam optimizer:

Initialize Adam optimizer parameters

eps = 1e-8
beta1 = 0.9
beta2 = 0.999
m_w = np.zeros(X_train.shape[1])
m_b = 0
v_w = np.zeros(X_train.shape[1])
v_b = 0

for i in range(num_iterations):
# Forward pass
z = np.dot(X_train, w) + b
y_pred = sigmoid(z)

# Compute the cost
cost = -np.mean(y_train * np.log(y_pred) + (1 - y_train) * np.log(1 - y_pred))

# Backward pass
dz = y_pred - y_train
dw = np.dot(X_train.T, dz) / X_train.shape[0]
db = np.mean(dz)

# Update the Adam optimizer parameters
m_w = beta1 * m_w + (1 - beta1) * dw
m_b = beta1 * m_b + (1 - beta1) * db
v_w = beta2 * v_w + (1 - beta2) * np.square(dw)
v_b = beta2 * v_b + (1 - beta2) * np.square(db)
m_w_hat = m_w / (1 - beta1**(i+1))
m_b_hat = m_b / (1 - beta1**(i+1))
v_w_hat = v_w / (1 - beta2**(i+1))
v_b_hat = v_b / (1 - beta2**(i+1))

# Update the parameters using Adam optimizer
w = w - learning_rate * m_w_hat / (np.sqrt(v_w_hat) + eps)
b = b - learning_rate * m_b_hat / (np.sqrt(v_b_hat) + eps)

Evaluate the model on the validation data

z_val = np.dot(X_val, w) + b
y_val_pred = sigmoid(z_val)
y_val_pred[y_val_pred >= 0.5] = 1
y_val_pred[y_val_pred < 0.5] = 0
accuracy = np.mean(y_val_pred == y_val)
print("Validation accuracy:", accuracy)

Propack currently prints to stdout. Can we trap these?

When evaluating Hessian in AMPL, I noticed that it returns a full sparse matrix, even though it is a symmetric matrix.
Is there a way to only have the lower triangular part by setting an option?

If not, I will do it in cython.

x = data[['variance', 'skewness']].values
y = data['class'].values
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(x, y)

X_train.shape, X_test.shape, y_train.shape, y_test.shape
x = (x - np.mean(x, axis=0)) / np.std(x, axis=0)
plt.scatter(x[:, 0], x[:, 1], c=y)
plt.show()

Perceptron algorithm

def perceptron(x, y, lr=0.1, n_iters=100):
w = np.zeros(x.shape[1])
b = 0

for _ in range(n_iters):
    for i in range(x.shape[0]):
        if y[i] * (np.dot(x[i], w) + b) <= 0:
            w += lr * y[i] * x[i]
            b += lr * y[i]

return w, b

perceptron_w, perceptron_b = perceptron(x, y)

Plot decision boundary for Perceptron

plt.scatter(x[:, 0], x[:, 1], c=y)
xlim = plt.gca().get_xlim()
ylim = plt.gca().get_ylim()
x_boundary = np.linspace(xlim[0], xlim[1])
y_boundary = -(perceptron_w[0] / perceptron_w[1]) * x_boundary - (perceptron_b / perceptron_w[1])
plt.plot(x_boundary, y_boundary, color='black')
plt.xlim(xlim)
plt.ylim(ylim)
plt.show()

Adaline algorithm

def adaline(x, y, lr=0.1, n_iters=100):
w = np.zeros(x.shape[1])
b = 0

for i in range(n_iters):
    output = np.dot(x, w) + b
    errors = y - output
    w += lr * np.dot(x.T, errors)
    b += lr * errors.sum()

return w, b

adaline_w, adaline_b = adaline(x, y)

Plot decision boundary for Adaline

plt.scatter(x[:, 0], x[:, 1], c=y)
xlim = plt.gca().get_xlim()
ylim = plt.gca().get_ylim()
x_boundary = np.linspace(xlim[0], xlim[1])
y_boundary = -(adaline_w[0] / adaline_w[1]) * x_boundary - (adaline_b / adaline_w[1])
plt.plot(x_boundary, y_boundary, color='black')
plt.xlim(xlim)
plt.ylim(ylim)
plt.show()

I'm planning to move lbfgs.py to PyKrylov. In Julia, the LBFGS operator is part of LinearOperators.jl and in Matlab, it's part of Spot. Does that sound reasonable? @syarra

import numpy as np
import matplotlib.pyplot as plt
from mlxtend.data import mnist_data

Load the MNIST dataset

X, y = mnist_data()

Shuffle the data

np.random.seed(42)
indices = np.random.permutation(len(X))
X = X[indices]
y = y[indices]

Standardize your dataset

X_mean = np.mean(X, axis=0)
X_std = np.std(X, axis=0)
X_std[X_std == 0] = 1e-8
X = (X - X_mean) / X_std

Divide data into training and test sets

train_ratio = 0.8 # 80% for training, 20% for testing
split_index = int(len(X) * train_ratio)
X_train = X[:split_index]
y_train = y[:split_index]
X_test = X[split_index:]
y_test = y[split_index:]

Apply one-hot encoding to the labels

num_classes = 10 # Number of classes in MNIST dataset
y_train_onehot = np.zeros((len(y_train), num_classes))
y_test_onehot = np.zeros((len(y_test), num_classes))
for i in range(len(y_train)):
y_train_onehot[i, y_train[i]] = 1
for i in range(len(y_test)):
y_test_onehot[i, y_test[i]] = 1

def sigmoid(z):
return 1 / (1 + np.exp(-z))

def mse_loss(y_true, y_pred):
return np.mean((y_true - y_pred) ** 2)

def calculate_accuracy(X, y, weights, biases):
outputs = forward_pass(X, weights, biases)
predicted_labels = np.argmax(outputs[-1], axis=1)
true_labels = np.argmax(y, axis=1)
accuracy = np.mean(predicted_labels == true_labels)
return accuracy

#------------------------------------------------------------------------------------------------------------

def initialize_parameters(input_size, hidden_layer_sizes, output_size):
# Combine input size, hidden layer sizes, and output size into a list
sizes = [input_size] + hidden_layer_sizes + [output_size]
weights = [] # Initialize an empty list to store the weights
biases = [] # Initialize an empty list to store the biases
for i in range(len(sizes) - 1): # Iterate over the sizes list, excluding the last element(output layer)
# Initialize random weights for the current layer
w = np.random.randn(sizes[i], sizes[i + 1]) # w [inputs from the previous layer * the corresponding weights.]
b = np.zeros((1, sizes[i + 1])) # Initialize biases as zeros for the current layer
weights.append(w) # Add the weights to the weights list
biases.append(b) # Add the biases to the biases list
return weights, biases # Return the weights and biases lists

def forward_pass(X, weights, biases):
# Store the input data as the first element of the outputs list
outputs = [X]
# Iterate over the weights list, representing each layer in the network
for i in range(len(weights)):
# Retrieve the outputs from the previous layer as the inputs to the current layer
inputs = outputs[-1]
# Compute the weighted sum of inputs and biases for the current layer
z = np.dot(inputs, weights[i]) + biases[i]
# Apply the sigmoid activation function to the weighted sum
a = sigmoid(z)
# Add the activations to the outputs list
outputs.append(a)
return outputs # Return the list of outputs for each layer

def backward_pass(X, y, weights, biases, outputs, learning_rate):
gradients = [] # Initialize an empty list to store the gradients of the weights.
num_samples = X.shape[0] # Get the number of samples in the input data.
d_output = (outputs[-1] - y) / num_samples # Calculate the gradient of the output layer activations.

# Compute the gradient of the weights connecting the last hidden layer to the output layer.
gradients.append(np.dot(outputs[-2].T, d_output))

biases_gradients = [np.mean(d_output, axis=0)]  # Compute the gradient of the biases for the output layer.

for i in range(len(weights) - 1, 0, -1):
    # Calculate the gradient of the hidden layer activations.
    d_hidden = np.dot(d_output, weights[i].T) * outputs[i] * (1 - outputs[i])

    # Compute the gradient of the weights connecting the previous layer to the current layer.
    gradients.append(np.dot(outputs[i - 1].T, d_hidden))
    # Compute the gradient of the biases for the current hidden layer.
    biases_gradients.append(np.mean(d_hidden, axis=0))

    d_output = d_hidden  # Update the gradient of the output layer activations.

gradients.reverse()  # Reverse the order of the gradients list.
biases_gradients.reverse()  # Reverse the order of the biases_gradients list.

for i in range(len(weights)):
    weights[i] -= learning_rate * gradients[i]  # Update the weights of each layer.
    biases[i] -= learning_rate * biases_gradients[i]  # Update the biases of each layer.

return weights, biases  # Return the updated weights and biases.

def train(X_train, y_train, X_test, y_test, num_of_layers, size_of_layers, learning_rate, num_epochs):
# Get the input and output sizes
input_size = X_train.shape[1]
output_size = y_train.shape[1]

# Initialize weights and biases
weights, biases = initialize_parameters(input_size, size_of_layers, output_size)

# Initialize lists to store losses and accuracies
train_losses = []
test_losses = []
accuracies = []

# Loop through each epoch
for epoch in range(num_epochs):
    # Forward pass on the training data
    train_outputs = forward_pass(X_train, weights, biases)

    # Calculate the training loss
    train_loss = mse_loss(y_train, train_outputs[-1])

    # Perform backward pass to update weights and biases
    weights, biases = backward_pass(X_train, y_train, weights, biases, train_outputs, learning_rate)

    # Store the training loss
    train_losses.append(train_loss)

    # Forward pass on the test data
    test_outputs = forward_pass(X_test, weights, biases)

    # Calculate the test loss
    test_loss = mse_loss(y_test, test_outputs[-1])

    # Store the test loss
    test_losses.append(test_loss)

    # Calculate accuracy on the test data
    accuracy = calculate_accuracy(X_test, y_test, weights, biases)

    # Store the accuracy
    accuracies.append(accuracy)

    # Print the progress every 500 epochs
    if (epoch + 1) % 500 == 0:
        print(
            f"Epoch {epoch + 1}/{num_epochs}, Train Loss: {train_loss:.4f}, Test Loss: {test_loss:.4f}, Accuracy: {accuracy:.4f}")

# Return the final weights, biases, losses, and accuracies
return weights, biases, train_losses, test_losses, accuracies

#-------------------------------------------------------------------------------------------------------------------

Neural Network function

def NN(X_train, y_train, X_test, y_test, num_of_layers, size_of_layers):
# Set the learning rate
learning_rate = 0.6

# Set the number of epochs
num_epochs = 1000

# Call the train function to train the neural network
weights, biases, train_losses, test_losses, accuracies = train(X_train, y_train, X_test, y_test, num_of_layers,
                                                               size_of_layers, learning_rate, num_epochs)

# Return the weights, biases, losses, and accuracies
return weights, biases, train_losses, test_losses, accuracies

Example usage

num_of_layers = 2
size_of_layers = [500,
500] # is a list that determines the number of neurons in each hidden layer of the neural network
weights, biases, train_losses, test_losses, accuracies = NN(X_train, y_train_onehot, X_test, y_test_onehot,
num_of_layers, size_of_layers)

Test with different architectures

architectures = [
(2, [100, 10]), # Architecture 1: 2 layers => 1 hidden layer (100 neurons) and 1 output layer (10 neurons)
(3, [50, 100, 10]), # Architecture 2: 3 layers => 2 hidden layers (100 and 50 neurons) and 1 output layer (10 neurons)
(3, [100, 50, 10]) # Architecture 3: 3 layers => 2 hidden layers (50 and 100 neurons) and 1 output layer (10 neurons)
]

for i, (num_of_layers, size_of_layers) in enumerate(architectures):
print(f"\nTesting Architecture {i+1}:")
weights, biases, train_losses, test_losses, accuracies = NN(X_train, y_train_onehot, X_test, y_test_onehot,
num_of_layers, size_of_layers)
final_accuracy = accuracies[-1]
print(f"Final Accuracy for Architecture {i+1}: {final_accuracy * 100}%")

#-------------------------------------------------------------------------------------------------------------------

Plot the training and test losses

plt.plot(train_losses, label='Training Loss')
plt.plot(test_losses, label='Test Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Training and Test Loss')
plt.legend()
plt.show()

final_accuracy = accuracies[-1]
print("Final Accuracy:", final_accuracy * 100)