The following example shows three different modes of applying archetypal analysis
to do dictionary learning. The input is 8 by 8 image patches extracted from one image.
The input patches are stored as columns of X. The learned archetypes are stored as columns of Z.

Load and Preprocess

Here we try to load the image 'lena.png' and to extract patches from it.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
import spams
import Image
import numpy as np

img_file = 'lena.png'
try:
    img = Image.open(img_file)
except Exception as e:
    print "Cannot load image %s (%s) : skipping test" %(img_file,e)
    raise
I = np.array(img) / 255.
if I.ndim == 3:
    A = np.asfortranarray(I.reshape((I.shape[0],I.shape[1] * I.shape[2])),dtype = np.float64)
    rgb = True
else:
    A = np.asfortranarray(I,dtype = np.float64)
    rgb = False

m = 8;n = 8;
X = spams.im2col_sliding(A,m,n,rgb)

The patches then are centered and normalized.

1
2
X = X - np.tile(np.mean(X,0),(X.shape[0],1))
X = np.asfortranarray(X / np.tile(np.sqrt((X * X).sum(axis=0)),(X.shape[0],1)),dtype = np.float64)

Set the parameters

Only K need to be set in the simplest setting. Other parameters are set as default if we don't specify them.

1
2
3
4
5
6
7
8
9
K = 64 # learns a dictionary with 64 elements
robust = False # use robust archetypal analysis or not, default parameter(True)
epsilon = 1e-3 # width in Huber loss, default parameter(1e-3)
computeXtX = True # memorize the product XtX or not default parameter(True)
stepsFISTA = 0 # 3 alternations by FISTA, default parameter(3)
# a for loop in FISTA is used, we stop at 50 iterations
# remember that we are not guarantee to descent in FISTA step if 50 is too small
stepsAS = 10 # 7 alternations by activeSet, default parameter(50)
randominit = True # random initilazation, default parameter(True)

First experiment

The first experiment show how to use archetypalAnalysis to learn p archetypes from input X.
It will run few steps using FISTA and then using active-set method.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# The mode with X and p known
# learn archetypes using activeSet method for each convex sub-problem
(Z,A,B) = spams.archetypalAnalysis(np.asfortranarray(X[:, :10000]), returnAB= True, p = K, \
    robust = robust, epsilon = epsilon, computeXtX = computeXtX,  stepsFISTA = stepsFISTA , stepsAS = stepsAS, numThreads = -1)

print 'Evaluating cost function...'
alpha = spams.decompSimplex(np.asfortranarray(X[:, :10000]),Z = Z, computeXtX = True, numThreads = -1)
xd = X[:,:10000] - Z * alpha
R = np.sum(xd*xd)
print "objective function: %f" %R

Second experiment

The second experiment is similar to the first one, except that we could give an initial value of the dictionary Z0.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# The mode with X and Z0 known
# learn archetypes using activeSet method for each convex sub-problem
Z2 = spams.archetypalAnalysis(np.asfortranarray(X[:, :10000]), Z0 = Z, robust = robust, \
    epsilon = epsilon, computeXtX = computeXtX , stepsFISTA = stepsFISTA,stepsAS = stepsAS, numThreads = -1)

print 'Evaluating cost function...'
alpha = spams.decompSimplex(np.asfortranarray(X[:, :10000]),Z = np.asfortranarray(Z2), computeXtX = True, numThreads = -1)
xd = X[:,:10000] - Z2 * alpha
R = np.sum(xd*xd)
print "objective function: %f" %R

Third experiment

The third experiment show how to use robust archetypal analysis.

1
2
(Z3,A3,B3) = spams.archetypalAnalysis(np.asfortranarray(X[:, :10000]), returnAB= True, p = K, \
    robust = True, epsilon = epsilon, computeXtX = computeXtX,  stepsFISTA = stepsFISTA , stepsAS = stepsAS, numThreads = -1)