The Data¶
In [6]:
mydata = demo_1_data()
Lloyd's Algorithm¶
Very Sensitive to initialization values
Converges, but no guarantees (esp in case of bad initializations)
No Guarantee about the results either
Really Fast (if no catastrophy)
In [6]:
demo_1_lloyd(mydata, 3)
Different Similarity Measures¶
The data does seem to have Gaussian Blobs
- The problem with the data is different
Different similarity metric will probably not give different results
- Except in case of relatively bad similarity metrics
In [50]:
demo_1_lloyd2(mydata, 3, ['euclidean', 'cityblock', 'mahalanobis'], 800)
demo_1_lloyd2(mydata, 3, ['euclidean', 'cityblock', 'mahalanobis'], 10)
Hartigan's Algorithm¶
Converges quickly
Relatively Robust
Still sensitive to initialization of classes
In [10]:
demo_1_hart(mydata, 3, [12, 42, 999]) # These took some time to choose
Smarter Way?¶
We couldn't figure out any smarter way, rather than :
only recalculate objective function for the current class
- Not reliable, esp when the data is disproportionate among classes
calculate the change in objective function incrementally
Halved the number of data points for which the distance is calculated, compared to Naive
- Does not fully utilize the potential, eg vectorization
m[ki] = (n[ki] / (n[ki] - 1)) * (m[ki] - x / n[ki])
normxkk = norm(x - m, axis=1)
eki = np.sum(norm(Xki - m[ki])) - normxkk[ki]
ediffki = eki - e[ki]
ediff = ediffki + normxkk
kw = np.argmin(ediff)
MacQueen's Algorithm¶
Convenient for streams
Sensitive to order of data
- Essentially, still sensitive to initialization
In [51]:
demo_1_macqueen(mydata, 3, [12, 556])
Run times!¶
In [18]:
np.random.seed(9000)
init_c = np.copy(mydata[np.random.choice(np.arange(mydata.shape[0]), size=3)])
demo_1_t()
Syllabus¶
Apply K-means
Apply Spectral Clustering
Apply Spectral Clustering using Andrew Ng's Alorithm
In [20]:
my_data = demo_2_1()
Apply K-means on Data¶
In [22]:
demo_2_2(my_data)
Apply Spectral Clustering on Data¶
- Get a good result at beta = 11
- By observation, we see that some edges points would be mis-judged as beta grows from 1 to 15
- The Upper half contains 100 points, and so is the lower half.
Play around the number of halfs¶
- See how the number of halfs changes
In [24]:
demo_2_4(my_data)
Exam the Laplacian Matrix¶
- S = exp(-beta* |x_i-x_j|^2) which is indepandent on the data order
- D = Sum(j to n)(Sij) if i=j which is depandent on the data order
- L = D - S
Shuffle the data order to see the result¶
- we would have a differnet beta or even unable to get one sometime. Sometimes we got a lot
- But we see that the upper half gathered close to y=0 line, while the lower half spread around.
In [26]:
demo_2_5(my_data)
Spectral Clustering using Andrew Ng's Algorithm¶
- On Spectral Clustering: Analysis and an algorithm by Andre Ng, etal.
- S = exp(-|x_i-x_j|^2 / 2*(sigma^2))
- D = diagonal matrix whose (i, i)-element is the sum of A's i-th row
- L = D^(-1/2)AD^(-1/2)
- FInd the k largest eigenvectors of L and normalized them into matrix X
- Treating X as new set of data, apply k-means clustering
- Based on the clustering result to cluster the original data
In [27]:
scan_rex.demo1(my_data,2, 0.05, 1, 0.1) #(my_data,k, sigma, start, end, step):
In [9]:
scan_rex.demo1(mydata.T, 3, 0.1, 1, 0.1) #(my_data,k, sigma, start, end, step):
Shuffle the data order to see the result¶
- We would have same L on differnet odering of data
- As prediction, the resulting sigma would not change
In [28]:
idx = np.arange(0,200,1)
np.random.shuffle(idx)
scan_rex.demo1(my_data[:,idx],2, 0.05, 1, 0.1) #(my_data,k, sigma, start, end, step):
In [31]:
res_pca, res_lda = demo_3(data_X, data_y, 2)
In [32]:
evals = res_lda[-1][0]
plt.plot(evals)
Out[32]:
In [33]:
res_pca, res_lda = demo_3(data_X, data_y, 2, use_abs_evals=True)
plt.figure()
plt.plot(res_lda[-1][0])
Out[33]:
Let's see for 3D¶
In [34]:
res_pca, res_lda = demo_3(data_X, data_y, 3, use_abs_evals=True)
Observations¶
PCA finds the axes with maximum variance for the whole data
unsupervised algorithm
performs better for fewer samples per class
LDA finds the axes for best separation between classes
- supervised algorithm
Both still have the underlying assumption of data having Gaussian Distribution
We all read the paper¶
In general, we expected:
method_1- to give the best 10th degree polygon,
- be fast,
- not face issues
method_2- to give okay result (conditionally)
- be slow,
- to face numerical issues
method_3- to give a good 10th degree polygon,
- be fast,
- not face issues
method_4- to give results a bit better than
method_2 - be the slowest (added transform step)
- not face issues in this particular case (may be)
- to give results a bit better than
Let's test the expectations!¶
In [36]:
def trsf(x, div=100):
return x / div
n = 10
def method_1(hgt=hgt, wgt=wgt, n=n, x=x):
# regression using ployfit
c = poly.polyfit(hgt, wgt, n)
y = poly.polyval(x, c)
return c, y
def method_2(hgt=hgt, wgt=wgt, n=n, x=x):
# regression using the Vandermonde matrix and pinv
X = poly.polyvander(hgt, n)
c = np.dot(la.pinv(X), wgt)
y = np.dot(poly.polyvander(x,n), c)
return c, y
def method_3(hgt=hgt, wgt=wgt, n=n, x=x):
# regression using the Vandermonde matrix and lstsq
X = poly.polyvander(hgt, n)
c = la.lstsq(X, wgt)[0]
y = np.dot(poly.polyvander(x,n), c)
return c, y
def method_4(hgt=hgt, wgt=wgt, n=n, x=x, div=100):
# regression on transformed data using the Vandermonde
# matrix and pinv
X = poly.polyvander(trsf(hgt, div=div), n)
c = np.dot(la.pinv(X), wgt)
y = np.dot(poly.polyvander(trsf(x, div=div),n), c)
return c, y
demo_4_plot([m()[1] for m in [method_1, method_2, method_3, method_4]],
[m for m in ["method_1", "method_2", "method_3", "method_4"]])
What happend?¶
Warning raised by polyfit
- abdullah: it is actually helpful
method_4gives almost as 'good' a result asmethod_1- can: the
trsffunction is scaling the values by 1/100 - cifong: problems with raising $10^2$ values to powers of 10 might be getting resolved
- abdullah: actually, no. Look at the coefficients. All are of order $10^{10}$
- can: the
In [37]:
coeffs = np.array([m()[0] for m in [method_1, method_2, method_3, method_4]])
print(coeffs.T)
- umut: Let's Time them
In [38]:
%timeit method_1()
%timeit method_2()
%timeit method_3()
%timeit method_4()
- abdullah: the dataset is too small. Run again with 1000 copies
- it is still in milliseconds
- can: don't go above 100,000
In [39]:
hgt_11 = np.tile(hgt, (1000))
wgt_11 = np.tile(wgt, (1000))
print(hgt_11.shape[0])
%timeit method_1(hgt_11, wgt_11)
%timeit method_2(hgt_11, wgt_11)
%timeit method_3(hgt_11, wgt_11)
%timeit method_4(hgt_11, wgt_11)
- cifong: Let me dive into the transformation
How do results change for method_4 for different divisors¶
In [40]:
ys = []
for d in [1, 10, 50, 100]:
ys.append(method_4(div=d)[1])
demo_4_plot(ys, [1, 10, 50, 100])
ys = []
for d in [100, 500, 1000, 1100]:
ys.append(method_4(div=d)[1])
demo_4_plot(ys, [100, 500, 1000, 1100])
The fitting dimention went up as the divisor goes up,¶
- The value in given data is 3 digit decimal,
- divisor = 100 make the value to 1 digit
Floating point issue in Python¶
- By Python Documentation
- 0.1(decimal) is presented as
- 0.00011001100110011001100110011001100110011001100110011010(Binary)
- Print (0.1)
- 0.1000000000000000055511151231257827021181583404541015625
- https://docs.python.org/2/tutorial/floatingpoint.html
In [41]:
print(np.finfo(np.float32))
print(np.finfo(np.float64))
Looking in to numpy.linalg.pinv¶
numpy.linalg.pinv(a, rcond=1e-15)rcondis the precision of float 64, whichlinalgsupports- we found this code inside
In [ ]:
a, wrap = _makearray(a)
_assertNoEmpty2d(a)
a = a.conjugate()
u, s, vt = svd(a, 0)
m = u.shape[0]
n = vt.shape[1]
cutoff = rcond*maximum.reduce(s)
for i in range(min(n, m)):
if s[i] > cutoff: #Suspicious part, set value to zero
s[i] = 1./s[i]
else:
s[i] = 0.;
res = dot(transpose(vt), multiply(s[:, newaxis], transpose(u)))
return wrap(res)
Test this part¶
- See how many cutsin different divisor
In [44]:
count = []
cutoffValue = []
inteNI = np.arange(0.001,1,0.001)
intePI = np.arange(1,1000,1)
divisor = np.append(inteNI, intePI)
for i in divisor:
X = poly.polyvander(trsf(hgt,i), n)
cutCount, cutoff = countPinvCutoff(X)
list.append(count, cutCount)
list.append(cutoffValue, cutoff)
fig = plt.figure(figsize=(18, 6))
plt.subplot(121)
plt.scatter(divisor, cutoffValue)
plt.subplot(122)
plt.scatter(divisor, count)
Out[44]:
Discussion¶
- 1.Ther trend mathes the degeneration of fitting function that
- the lower the cut off, the higher remaining dimension we have
- 2.We would have at least 2 cutoffs in all divisor
In [45]:
count = []
for i in np.arange(0.0001,1,0.0001):
X = poly.polyvander(trsf(hgt,i), n)
count.append(countPinvCutoff(X)[0])
plt.scatter(np.arange(0.0001,1,0.0001), count)
Out[45]:
Plot some divisors on 3 variants¶
- Compare with the tolerence = epslon max(m,n) maximum(s)
- linalg.inv cutoff = epslon * maximum(s)
- No cut off
In [47]:
index = [0.01,0.05,0.1,0.5, 1,5, 10, 50, 100, 500, 1000]
for i in index:
demoPic(i)
