MNISTでソフトマックス回帰したいがうまくいかない

Question

前提・実現したいこと

MNISTで60000個の訓練データと10000個のテストデータを用いて、
ソフトマックス回帰を実装したいです。
ただしkerasはデータの読み込み以外に使わず、
ネットワーク自体は書きたいです(大学の課題なので)

以下のコードを実行しても、エラーは出ませんが計算が2時間以上終わりません。
以下のコードにまずい点はあるでしょうか？

該当のソースコード

import numpy as np
from keras.datasets import mnist # データの読み込み
from keras.utils import np_utils # 1-of-K符号化法
from scipy.optimize import minimize # 検算用

(x_train, y_train), (x_test, y_test) = mnist.load_data()

x_train = x_train.reshape(60000, 784) 
x_train = x_train.astype('float32')
x_train = x_train / 255
num_classes = 10
y_train = np_utils.to_categorical(y_train, num_classes)

x_test = x_test.reshape(10000, 784)
x_test = x_test.astype('float32')
x_test = x_test / 255
y_test = np_utils.to_categorical(y_test, num_classes)


# softmax関数
def logistic(x, w):
    w = w.reshape((10, 785)) # (クラス数 = 10, 次元数 + 1 = 785)
    x = np.c_[x, np.ones((60000, 1))] # x_trainの最後列に1を追加する
    y = np.exp(np.dot(x, w.T)) # y = exp(線形和)

    wk = np.sum(y, axis = 1) 
    wk = y.T / wk
    y = wk.T
    return y

# 交差エントロピー誤差
def cee_logistic(w, x, t):
    y = logistic(x, w)
    y = np.log(y) # 先にlogを取っておく
    N = y.shape[0]

    cee = 0
    for n in range(N):
        cee = cee - np.dot(t[n, :], y[n, :])
    
    cee = cee / N
    return cee

# 交差エントロピー誤差の微分
def dcee_logistic(w, x, t):
    y = logistic(x, w)
    N, K = y.shape
    x = np.c_[x, np.ones((60000, 1))] # x_trainの最後列に1を追加する
    
    dcee = np.zeros((10, 785)) #交差エントロピー誤差を重みで微分した行列
    for k in range(K):
        wk = y[:, k] - t[:, k]
        wk = np.dot(x.T ,wk)
        dcee[k, :] = wk.T

    dcee = dcee / N
    return dcee.reshape(-1)

# 最適化
def fit_logistic(w_init, x, t):
    res = minimize(fun = cee_logistic, x0 = w_init, args = (x, t), jac = dcee_logistic, method = "CG")
    return res.x

W_init = np.zeros((10, 785))
W = fit_logistic(W_init, x_train, y_train)
print(np.round(W.reshape(10, 785), 2))