Trial 2: classification with learned graph filters

We want to classify data by first extracting meaningful features from learned filters.


In [ ]:
import time
import numpy as np
import scipy.sparse, scipy.sparse.linalg, scipy.spatial.distance
from sklearn import datasets, linear_model
import matplotlib.pyplot as plt
%matplotlib inline

import os
import sys
sys.path.append('..')
from lib import graph

Parameters

Dataset

  • Two digits version of MNIST with N samples of each class.
  • Distinguishing 4 from 9 is the hardest.

In [ ]:
def mnist(a, b, N):
    """Prepare data for binary classification of MNIST."""
    folder = os.path.join('..', 'data')
    mnist = datasets.fetch_mldata('MNIST original', data_home=folder)

    assert N < min(sum(mnist.target==a), sum(mnist.target==b))
    M = mnist.data.shape[1]
    
    X = np.empty((M, 2, N))
    X[:,0,:] = mnist.data[mnist.target==a,:][:N,:].T
    X[:,1,:] = mnist.data[mnist.target==b,:][:N,:].T
    
    y = np.empty((2, N))
    y[0,:] = -1
    y[1,:] = +1

    X.shape = M, 2*N
    y.shape = 2*N, 1
    return X, y

X, y = mnist(4, 9, 1000)

print('Dimensionality: N={} samples, M={} features'.format(X.shape[1], X.shape[0]))

X -= 127.5
print('X in [{}, {}]'.format(np.min(X), np.max(X)))

def plot_digit(nn):
    M, N = X.shape
    m = int(np.sqrt(M))
    fig, axes = plt.subplots(1,len(nn), figsize=(15,5))
    for i, n in enumerate(nn):
        n = int(n)
        img = X[:,n]
        axes[i].imshow(img.reshape((m,m)))
        axes[i].set_title('Label: y = {:.0f}'.format(y[n,0]))

plot_digit([0, 1, 1e2, 1e2+1, 1e3, 1e3+1])

Regularized least-square

Reference: sklearn ridge regression

  • With regularized data, the objective is the same with or without bias.

In [ ]:
def test_sklearn(tauR):
    
    def L(w, b=0):
        return np.linalg.norm(X.T @ w + b - y)**2 + tauR * np.linalg.norm(w)**2

    def dL(w):
        return 2 * X @ (X.T @ w - y) + 2 * tauR * w

    clf = linear_model.Ridge(alpha=tauR, fit_intercept=False)
    clf.fit(X.T, y)
    w = clf.coef_.T

    print('L = {}'.format(L(w, clf.intercept_)))
    print('|dLw| = {}'.format(np.linalg.norm(dL(w))))

    # Normalized data: intercept should be small.
    print('bias: {}'.format(abs(np.mean(y - X.T @ w))))

test_sklearn(1e-3)

Linear classifier


In [ ]:
def test_optim(clf, X, y, ax=None):
    """Test optimization on full dataset."""
    tstart = time.process_time()
    ret = clf.fit(X, y)
    print('Processing time: {}'.format(time.process_time()-tstart))
    print('L = {}'.format(clf.L(*ret, y)))
    if hasattr(clf, 'dLc'):
        print('|dLc| = {}'.format(np.linalg.norm(clf.dLc(*ret, y))))
    if hasattr(clf, 'dLw'):
        print('|dLw| = {}'.format(np.linalg.norm(clf.dLw(*ret, y))))
    if hasattr(clf, 'loss'):
        if not ax:
            fig = plt.figure()
            ax = fig.add_subplot(111)
        ax.semilogy(clf.loss)
        ax.set_title('Convergence')
        ax.set_xlabel('Iteration number')
        ax.set_ylabel('Loss')
    if hasattr(clf, 'Lsplit'):
        print('Lsplit = {}'.format(clf.Lsplit(*ret, y)))
        print('|dLz| = {}'.format(np.linalg.norm(clf.dLz(*ret, y))))
        ax.semilogy(clf.loss_split)

In [ ]:
class rls:
    
    def __init__(s, tauR, algo='solve'):
        s.tauR = tauR
        if algo is 'solve':
            s.fit = s.solve
        elif algo is 'inv':
            s.fit = s.inv

    def L(s, X, y):
        return np.linalg.norm(X.T @ s.w - y)**2 + s.tauR * np.linalg.norm(s.w)**2

    def dLw(s, X, y):
        return 2 * X @ (X.T @ s.w - y) + 2 * s.tauR * s.w
    
    def inv(s, X, y):
        s.w = np.linalg.inv(X @ X.T + s.tauR * np.identity(X.shape[0])) @ X @ y
        return (X,)
    
    def solve(s, X, y):
        s.w = np.linalg.solve(X @ X.T + s.tauR * np.identity(X.shape[0]), X @ y)
        return (X,)
    
    def predict(s, X):
        return X.T @ s.w

test_optim(rls(1e-3, 'solve'), X, y)
test_optim(rls(1e-3, 'inv'), X, y)

Feature graph


In [ ]:
t_start = time.process_time()
z = graph.grid(int(np.sqrt(X.shape[0])))
dist, idx = graph.distance_sklearn_metrics(z, k=4)
A = graph.adjacency(dist, idx)
L = graph.laplacian(A, True)
lmax = graph.lmax(L)
print('Execution time: {:.2f}s'.format(time.process_time() - t_start))

Lanczos basis


In [ ]:
def lanczos(L, X, K):
    M, N = X.shape
    a = np.empty((K, N))
    b = np.zeros((K, N))
    V = np.empty((K, M, N))
    V[0,...] = X / np.linalg.norm(X, axis=0)
    for k in range(K-1):
        W = L.dot(V[k,...])
        a[k,:] = np.sum(W * V[k,...], axis=0)
        W = W - a[k,:] * V[k,...] - (b[k,:] * V[k-1,...] if k>0 else 0)
        b[k+1,:] = np.linalg.norm(W, axis=0)
        V[k+1,...] = W / b[k+1,:]
    a[K-1,:] = np.sum(L.dot(V[K-1,...]) * V[K-1,...], axis=0)
    return V, a, b

def lanczos_H_diag(a, b):
    K, N = a.shape
    H = np.zeros((K*K, N))
    H[:K**2:K+1, :] = a
    H[1:(K-1)*K:K+1, :] = b[1:,:]
    H.shape = (K, K, N)
    Q = np.linalg.eigh(H.T, UPLO='L')[1]
    Q = np.swapaxes(Q,1,2).T
    return Q

def lanczos_basis_eval(L, X, K):
    V, a, b = lanczos(L, X, K)
    Q = lanczos_H_diag(a, b)
    M, N = X.shape
    Xt = np.empty((K, M, N))
    for n in range(N):
        Xt[...,n] = Q[...,n].T @ V[...,n]
    Xt *= Q[0,:,np.newaxis,:]
    Xt *= np.linalg.norm(X, axis=0)
    return Xt, Q[0,...]

Tests

  • Memory arrangement for fastest computations: largest dimensions on the outside, i.e. fastest varying indices.
  • The einsum seems to be efficient for three operands.

In [ ]:
def test():
    """Test the speed of filtering and weighting."""
    
    def mult(impl=3):
        if impl is 0:
            Xb = Xt.view()
            Xb.shape = (K, M*N)
            XCb = Xb.T @ C  # in MN x F
            XCb = XCb.T.reshape((F*M, N))
            return (XCb.T @ w).squeeze()
        elif impl is 1:
            tmp = np.tensordot(Xt, C, (0,0))
            return np.tensordot(tmp, W, ((0,2),(1,0)))
        elif impl is 2:
            tmp = np.tensordot(Xt, C, (0,0))
            return np.einsum('ijk,ki->j', tmp, W)
        elif impl is 3:
            return np.einsum('kmn,fm,kf->n', Xt, W, C)
    
    C = np.random.normal(0,1,(K,F))
    W = np.random.normal(0,1,(F,M))
    w = W.reshape((F*M, 1))
    a = mult(impl=0)
    for impl in range(4):
        tstart = time.process_time()
        for k in range(1000):
            b = mult(impl)
        print('Execution time (impl={}): {}'.format(impl, time.process_time() - tstart))
        np.testing.assert_allclose(a, b)
#test()

GFL classification without weights

  • The matrix is singular thus not invertible.

In [ ]:
class gflc_noweights:

    def __init__(s, F, K, niter, algo='direct'):
        """Model hyper-parameters"""
        s.F = F
        s.K = K
        s.niter = niter
        if algo is 'direct':
            s.fit = s.direct
        elif algo is 'sgd':
            s.fit = s.sgd
    
    def L(s, Xt, y):
        #tmp = np.einsum('kmn,kf,fm->n', Xt, s.C, np.ones((s.F,M))) - y.squeeze()
        #tmp = np.einsum('kmn,kf->mnf', Xt, s.C).sum((0,2)) - y.squeeze()
        #tmp = (C.T @ Xt.reshape((K,M*N))).reshape((F,M,N)).sum((0,2)) - y.squeeze()
        tmp = np.tensordot(s.C, Xt, (0,0)).sum((0,1)) - y.squeeze()
        return np.linalg.norm(tmp)**2

    def dLc(s, Xt, y):
        tmp = np.tensordot(s.C, Xt, (0,0)).sum(axis=(0,1)) - y.squeeze()
        return np.dot(Xt, tmp).sum(1)[:,np.newaxis].repeat(s.F,1)
        #return np.einsum('kmn,n->km', Xt, tmp).sum(1)[:,np.newaxis].repeat(s.F,1)

    def sgd(s, X, y):
        Xt, q = lanczos_basis_eval(L, X, s.K)
        s.C = np.random.normal(0, 1, (s.K, s.F))
        s.loss = [s.L(Xt, y)]
        for t in range(s.niter):
            s.C -= 1e-13 * s.dLc(Xt, y)
            s.loss.append(s.L(Xt, y))
        return (Xt,)
    
    def direct(s, X, y):
        M, N = X.shape
        Xt, q = lanczos_basis_eval(L, X, s.K)
        s.C = np.random.normal(0, 1, (s.K, s.F))
        W = np.ones((s.F, M))
        c = s.C.reshape((s.K*s.F, 1))
        s.loss = [s.L(Xt, y)]
        Xw = np.einsum('kmn,fm->kfn', Xt, W)
        #Xw = np.tensordot(Xt, W, (1,1))
        Xw.shape = (s.K*s.F, N)
        #np.linalg.inv(Xw @ Xw.T)
        c[:] = np.linalg.solve(Xw @ Xw.T, Xw @ y)
        s.loss.append(s.L(Xt, y))
        return (Xt,)

#test_optim(gflc_noweights(1, 4, 100, 'sgd'), X, y)
#test_optim(gflc_noweights(1, 4, 0, 'direct'), X, y)

GFL classification with weights


In [ ]:
class gflc_weights():

    def __init__(s, F, K, tauR, niter, algo='direct'):
        """Model hyper-parameters"""
        s.F = F
        s.K = K
        s.tauR = tauR
        s.niter = niter
        if algo is 'direct':
            s.fit = s.direct
        elif algo is 'sgd':
            s.fit = s.sgd

    def L(s, Xt, y):
        tmp = np.einsum('kmn,kf,fm->n', Xt, s.C, s.W) - y.squeeze()
        return np.linalg.norm(tmp)**2 + s.tauR * np.linalg.norm(s.W)**2

    def dLw(s, Xt, y):
        tmp = np.einsum('kmn,kf,fm->n', Xt, s.C, s.W) - y.squeeze()
        return 2 * np.einsum('kmn,kf,n->fm', Xt, s.C, tmp) + 2 * s.tauR * s.W

    def dLc(s, Xt, y):
        tmp = np.einsum('kmn,kf,fm->n', Xt, s.C, s.W) - y.squeeze()
        return 2 * np.einsum('kmn,n,fm->kf', Xt, tmp, s.W)

    def sgd(s, X, y):
        M, N = X.shape
        Xt, q = lanczos_basis_eval(L, X, s.K)
        s.C = np.random.normal(0, 1, (s.K, s.F))
        s.W = np.random.normal(0, 1, (s.F, M))

        s.loss = [s.L(Xt, y)]

        for t in range(s.niter):
            s.C -= 1e-12 * s.dLc(Xt, y)
            s.W -= 1e-12 * s.dLw(Xt, y)
            s.loss.append(s.L(Xt, y))
        
        return (Xt,)

    def direct(s, X, y):
        M, N = X.shape
        Xt, q = lanczos_basis_eval(L, X, s.K)
        s.C = np.random.normal(0, 1, (s.K, s.F))
        s.W = np.random.normal(0, 1, (s.F, M))
        #c = s.C.reshape((s.K*s.F, 1))
        #w = s.W.reshape((s.F*M, 1))
        c = s.C.view()
        c.shape = (s.K*s.F, 1)
        w = s.W.view()
        w.shape = (s.F*M, 1)

        s.loss = [s.L(Xt, y)]

        for t in range(s.niter):
            Xw = np.einsum('kmn,fm->kfn', Xt, s.W)
            #Xw = np.tensordot(Xt, s.W, (1,1))
            Xw.shape = (s.K*s.F, N)
            c[:] = np.linalg.solve(Xw @ Xw.T, Xw @ y)

            Z = np.einsum('kmn,kf->fmn', Xt, s.C)
            #Z = np.tensordot(Xt, s.C, (0,0))
            #Z = s.C.T @ Xt.reshape((K,M*N))
            Z.shape = (s.F*M, N)
            w[:] = np.linalg.solve(Z @ Z.T + s.tauR * np.identity(s.F*M), Z @ y)

            s.loss.append(s.L(Xt, y))
        
        return (Xt,)

    def predict(s, X):
        Xt, q = lanczos_basis_eval(L, X, s.K)
        return np.einsum('kmn,kf,fm->n', Xt, s.C, s.W)

#test_optim(gflc_weights(3, 4, 1e-3, 50, 'sgd'), X, y)
clf_weights = gflc_weights(F=3, K=50, tauR=1e4, niter=5, algo='direct')
test_optim(clf_weights, X, y)

GFL classification with splitting

Solvers

  • Closed-form solution.
  • Stochastic gradient descent.

In [ ]:
class gflc_split():

    def __init__(s, F, K, tauR, tauF, niter, algo='direct'):
        """Model hyper-parameters"""
        s.F = F
        s.K = K
        s.tauR = tauR
        s.tauF = tauF
        s.niter = niter
        if algo is 'direct':
            s.fit = s.direct
        elif algo is 'sgd':
            s.fit = s.sgd

    def L(s, Xt, XCb, Z, y):
        return np.linalg.norm(XCb.T @ s.w - y)**2 + s.tauR * np.linalg.norm(s.w)**2

    def Lsplit(s, Xt, XCb, Z, y):
        return np.linalg.norm(Z.T @ s.w - y)**2 + s.tauF * np.linalg.norm(XCb - Z)**2 + s.tauR * np.linalg.norm(s.w)**2

    def dLw(s, Xt, XCb, Z, y):
        return 2 * Z @ (Z.T @ s.w - y) + 2 * s.tauR * s.w

    def dLc(s, Xt, XCb, Z, y):
        Xb = Xt.reshape((s.K, -1)).T
        Zb = Z.reshape((s.F, -1)).T
        return 2 * s.tauF * Xb.T @ (Xb @ s.C - Zb)

    def dLz(s, Xt, XCb, Z, y):
        return 2 * s.w @ (s.w.T @ Z - y.T) + 2 * s.tauF * (Z - XCb)

    def lanczos_filter(s, Xt):
        M, N = Xt.shape[1:]
        Xb = Xt.reshape((s.K, M*N)).T
        #XCb = np.tensordot(Xb, C, (2,1))
        XCb = Xb @ s.C  # in MN x F
        XCb = XCb.T.reshape((s.F*M, N))  # Needs to copy data.
        return XCb

    def sgd(s, X, y):
        M, N = X.shape
        Xt, q = lanczos_basis_eval(L, X, s.K)
        s.C = np.zeros((s.K, s.F))
        s.w = np.zeros((s.F*M, 1))
        Z = np.random.normal(0, 1, (s.F*M, N))

        XCb = np.empty((s.F*M, N))

        s.loss = [s.L(Xt, XCb, Z, y)]
        s.loss_split = [s.Lsplit(Xt, XCb, Z, y)]

        for t in range(s.niter):
            s.C -= 1e-7 * s.dLc(Xt, XCb, Z, y)
            XCb[:] = s.lanczos_filter(Xt)
            Z -= 1e-4 * s.dLz(Xt, XCb, Z, y)
            s.w -= 1e-4 * s.dLw(Xt, XCb, Z, y)
            s.loss.append(s.L(Xt, XCb, Z, y))
            s.loss_split.append(s.Lsplit(Xt, XCb, Z, y))
        
        return Xt, XCb, Z

    def direct(s, X, y):
        M, N = X.shape
        Xt, q = lanczos_basis_eval(L, X, s.K)
        s.C = np.zeros((s.K, s.F))
        s.w = np.zeros((s.F*M, 1))
        Z = np.random.normal(0, 1, (s.F*M, N))

        XCb = np.empty((s.F*M, N))
        Xb = Xt.reshape((s.K, M*N)).T
        Zb = Z.reshape((s.F, M*N)).T

        s.loss = [s.L(Xt, XCb, Z, y)]
        s.loss_split = [s.Lsplit(Xt, XCb, Z, y)]

        for t in range(s.niter):

            s.C[:] = Xb.T @ Zb / np.sum((np.linalg.norm(X, axis=0) * q)**2, axis=1)[:,np.newaxis]
            XCb[:] = s.lanczos_filter(Xt)

            #Z[:] = np.linalg.inv(s.tauF * np.identity(s.F*M) + s.w @ s.w.T) @ (s.tauF * XCb + s.w @ y.T)
            Z[:] = np.linalg.solve(s.tauF * np.identity(s.F*M) + s.w @ s.w.T, s.tauF * XCb + s.w @ y.T)

            #s.w[:] = np.linalg.inv(Z @ Z.T + s.tauR * np.identity(s.F*M)) @ Z @ y
            s.w[:] = np.linalg.solve(Z @ Z.T + s.tauR * np.identity(s.F*M), Z @ y)

            s.loss.append(s.L(Xt, XCb, Z, y))
            s.loss_split.append(s.Lsplit(Xt, XCb, Z, y))
        
        return Xt, XCb, Z

    def predict(s, X):
        Xt, q = lanczos_basis_eval(L, X, s.K)
        XCb = s.lanczos_filter(Xt)
        return XCb.T @ s.w

#test_optim(gflc_split(3, 4, 1e-3, 1e-3, 50, 'sgd'), X, y)
clf_split = gflc_split(3, 4, 1e4, 1e-3, 8, 'direct')
test_optim(clf_split, X, y)

Filters visualization

Observations:

  • Filters learned with the splitting scheme have much smaller amplitudes.
  • Maybe the energy sometimes goes in W ?
  • Why are the filters so different ?

In [ ]:
lamb, U = graph.fourier(L)
print('Spectrum in [{:1.2e}, {:1.2e}]'.format(lamb[0], lamb[-1]))

In [ ]:
def plot_filters(C, spectrum=False):
    K, F = C.shape
    M, M = L.shape
    m = int(np.sqrt(M))
    X = np.zeros((M,1))
    X[int(m/2*(m+1))] = 1  # Kronecker
    Xt, q = lanczos_basis_eval(L, X, K)
    Z = np.einsum('kmn,kf->mnf', Xt, C)
    Xh = U.T @ X
    Zh = np.tensordot(U.T, Z, (1,0))
    
    pmin = int(m/2) - K
    pmax = int(m/2) + K + 1
    fig, axes = plt.subplots(2,int(np.ceil(F/2)), figsize=(15,5))
    for f in range(F):
        img = Z[:,0,f].reshape((m,m))[pmin:pmax,pmin:pmax]
        im = axes.flat[f].imshow(img, vmin=Z.min(), vmax=Z.max(), interpolation='none')
        axes.flat[f].set_title('Filter {}'.format(f))
    fig.subplots_adjust(right=0.8)
    cax = fig.add_axes([0.82, 0.16, 0.02, 0.7])
    fig.colorbar(im, cax=cax)
    
    if spectrum:
        ax = plt.figure(figsize=(15,5)).add_subplot(111)
        for f in range(F):
            ax.plot(lamb, Zh[...,f] / Xh, '.-', label='Filter {}'.format(f))
        ax.legend(loc='best')
        ax.set_title('Spectrum of learned filters')
        ax.set_xlabel('Frequency')
        ax.set_ylabel('Amplitude')
        ax.set_xlim(0, lmax)

plot_filters(clf_weights.C, True)
plot_filters(clf_split.C, True)

Extracted features


In [ ]:
def plot_features(C, x):
    K, F = C.shape
    m = int(np.sqrt(x.shape[0]))
    xt, q = lanczos_basis_eval(L, x, K)
    Z = np.einsum('kmn,kf->mnf', xt, C)
    
    fig, axes = plt.subplots(2,int(np.ceil(F/2)), figsize=(15,5))
    for f in range(F):
        img = Z[:,0,f].reshape((m,m))
        #im = axes.flat[f].imshow(img, vmin=Z.min(), vmax=Z.max(), interpolation='none')
        im = axes.flat[f].imshow(img, interpolation='none')
        axes.flat[f].set_title('Filter {}'.format(f))
    fig.subplots_adjust(right=0.8)
    cax = fig.add_axes([0.82, 0.16, 0.02, 0.7])
    fig.colorbar(im, cax=cax)

plot_features(clf_weights.C, X[:,[0]])
plot_features(clf_weights.C, X[:,[1000]])

Performance w.r.t. hyper-parameters

  • F plays a big role.
    • Both for performance and training time.
    • Larger values lead to over-fitting !
  • Order $K \in [3,5]$ seems sufficient.
  • $\tau_R$ does not have much influence.

In [ ]:
def scorer(clf, X, y):
    yest = clf.predict(X).round().squeeze()
    y = y.squeeze()
    yy = np.ones(len(y))
    yy[yest < 0] = -1
    nerrs = np.count_nonzero(y - yy)
    return 1 - nerrs / len(y)

In [ ]:
def perf(clf, nfolds=3):
    """Test training accuracy."""
    N = X.shape[1]
    inds = np.arange(N)
    np.random.shuffle(inds)
    inds.resize((nfolds, int(N/nfolds)))
    folds = np.arange(nfolds)
    test = inds[0,:]
    train = inds[folds != 0, :].reshape(-1)
    
    fig, axes = plt.subplots(1,3, figsize=(15,5))
    test_optim(clf, X[:,train], y[train], axes[2])
    
    axes[0].plot(train, clf.predict(X[:,train]), '.')
    axes[0].plot(train, y[train].squeeze(), '.')
    axes[0].set_ylim([-3,3])
    axes[0].set_title('Training set accuracy: {:.2f}'.format(scorer(clf, X[:,train], y[train])))
    axes[1].plot(test, clf.predict(X[:,test]), '.')
    axes[1].plot(test, y[test].squeeze(), '.')
    axes[1].set_ylim([-3,3])
    axes[1].set_title('Testing set accuracy: {:.2f}'.format(scorer(clf, X[:,test], y[test])))
    
    if hasattr(clf, 'C'):
        plot_filters(clf.C)

perf(rls(tauR=1e6))
for F in [1,3,5]:
    perf(gflc_weights(F=F, K=50, tauR=1e4, niter=5, algo='direct'))

#perf(rls(tauR=1e-3))
#for K in [2,3,5,7]:
#    perf(gflc_weights(F=3, K=K, tauR=1e-3, niter=5, algo='direct'))

#for tauR in [1e-3, 1e-1, 1e1]:
#    perf(rls(tauR=tauR))
#    perf(gflc_weights(F=3, K=3, tauR=tauR, niter=5, algo='direct'))

Classification

  • Greater is $F$, greater should $K$ be.

In [ ]:
def cross_validation(clf, nfolds, nvalidations):
    M, N = X.shape
    scores = np.empty((nvalidations, nfolds))
    for nval in range(nvalidations):
        inds = np.arange(N)
        np.random.shuffle(inds)
        inds.resize((nfolds, int(N/nfolds)))
        folds = np.arange(nfolds)
        for n in folds:
            test = inds[n,:]
            train = inds[folds != n, :].reshape(-1)
            clf.fit(X[:,train], y[train])
            scores[nval, n] = scorer(clf, X[:,test], y[test])
    return scores.mean()*100, scores.std()*100
    #print('Accuracy: {:.2f} +- {:.2f}'.format(scores.mean()*100, scores.std()*100))
    #print(scores)

In [ ]:
def test_classification(clf, params, param, values, nfolds=10, nvalidations=1):
    means = []
    stds = []
    fig, ax = plt.subplots(1,1, figsize=(15,5))
    for i,val in enumerate(values):
        params[param] = val
        mean, std = cross_validation(clf(**params), nfolds, nvalidations)
        means.append(mean)
        stds.append(std)
        ax.annotate('{:.2f} +- {:.2f}'.format(mean,std), xy=(i,mean), xytext=(10,10), textcoords='offset points')
    ax.errorbar(np.arange(len(values)), means, stds, fmt='.', markersize=10)
    ax.set_xlim(-.8, len(values)-.2)
    ax.set_xticks(np.arange(len(values)))
    ax.set_xticklabels(values)
    ax.set_xlabel(param)
    ax.set_ylim(50, 100)
    ax.set_ylabel('Accuracy')
    ax.set_title('Parameters: {}'.format(params))

In [ ]:
test_classification(rls, {}, 'tauR', [1e8,1e7,1e6,1e5,1e4,1e3,1e-5,1e-8], 10, 10)

In [ ]:
params = {'F':1, 'K':2, 'tauR':1e3, 'niter':5, 'algo':'direct'}
test_classification(gflc_weights, params, 'tauR', [1e8,1e6,1e5,1e4,1e3,1e2,1e-3,1e-8], 10, 10)

In [ ]:
params = {'F':2, 'K':10, 'tauR':1e4, 'niter':5, 'algo':'direct'}
test_classification(gflc_weights, params, 'F', [1,2,3,5])

In [ ]:
params = {'F':2, 'K':4, 'tauR':1e4, 'niter':5, 'algo':'direct'}
test_classification(gflc_weights, params, 'K', [2,3,4,5,8,10,20,30,50,70])

Sampled MNIST


In [ ]:
Xfull = X

In [ ]:
def sample(X, p, seed=None):
    M, N = X.shape
    z = graph.grid(int(np.sqrt(M)))
    
    # Select random pixels.
    np.random.seed(seed)
    mask = np.arange(M)
    np.random.shuffle(mask)
    mask = mask[:int(p*M)]
    
    return z[mask,:], X[mask,:]

X = Xfull
z, X = sample(X, .5)
dist, idx = graph.distance_sklearn_metrics(z, k=4)
A = graph.adjacency(dist, idx)
L = graph.laplacian(A)
lmax = graph.lmax(L)
lamb, U = graph.fourier(L)
print('Spectrum in [{:1.2e}, {:1.2e}]'.format(lamb[0], lamb[-1]))

print(L.shape)

def plot(n):
    M, N = X.shape
    m = int(np.sqrt(M))
    x = X[:,n]
    #print(x+127.5)
    plt.scatter(z[:,0], -z[:,1], s=20, c=x+127.5)
plot(10)

def plot_digit(nn):
    M, N = X.shape
    m = int(np.sqrt(M))
    fig, axes = plt.subplots(1,len(nn), figsize=(15,5))
    for i, n in enumerate(nn):
        n = int(n)
        img = X[:,n]
        axes[i].imshow(img.reshape((m,m)))
        axes[i].set_title('Label: y = {:.0f}'.format(y[n,0]))

#plot_digit([0, 1, 1e2, 1e2+1, 1e3, 1e3+1])

In [ ]:
#clf_weights = gflc_weights(F=3, K=4, tauR=1e-3, niter=5, algo='direct')
#test_optim(clf_weights, X, y)
#plot_filters(clf_weights.C, True)

In [ ]:
#test_classification(rls, {}, 'tauR', [1e1,1e0])
#params = {'F':2, 'K':5, 'tauR':1e-3, 'niter':5, 'algo':'direct'}
#test_classification(gflc_weights, params, 'F', [1,2,3])

In [ ]:
test_classification(rls, {}, 'tauR', [1e8,1e7,1e6,1e5,1e4,1e3,1e-5,1e-8], 10, 10)

In [ ]:
params = {'F':2, 'K':2, 'tauR':1e3, 'niter':5, 'algo':'direct'}
test_classification(gflc_weights, params, 'tauR', [1e8,1e5,1e4,1e3,1e2,1e1,1e-3,1e-8], 10, 1)

In [ ]:
params = {'F':2, 'K':10, 'tauR':1e5, 'niter':5, 'algo':'direct'}
test_classification(gflc_weights, params, 'F', [1,2,3,4,5,10])

In [ ]:
params = {'F':2, 'K':4, 'tauR':1e5, 'niter':5, 'algo':'direct'}
test_classification(gflc_weights, params, 'K', [2,3,4,5,6,7,8,10,20,30])