In [ ]:
function pca(x, n_comp, whitened)
-- Center data
local mean = torch.mean(x, 1)
local x_m = x - torch.ger(torch.ones(x:size(1)), mean:squeeze())
-- Calculate Covariance
local cov = x_m * x_m:t()
-- cov:div(x:size(1) - 1)
-- Get eigenvalues and eigenvectors
local ce, cv = torch.symeig(cov, 'V')
-- Sort eigenvalues
local ce, idx = torch.sort(ce, true)
-- Sort eigenvectors
cv = cv:index(2, idx:long())
-- Keep only the top
if n_comp and n_comp < cv:size(2) then
ce = ce:sub(1, n_comp)
cv = cv:sub(1, -1, 1, n_comp)
end
-- Check if whitened version
-- vectors are divided by the singular values to ensure uncorrelated outputs with unit component-wise variances.
if not whitened then
ce:add(1e-5):sqrt()
end
-- Get inverse
local inv_ce = ce:clone():pow(-1)
-- Make it a matrix with diagonal inv_ce
local inv_diag = torch.diag(inv_ce)
-- Transform to get U
local u = x_m:t() * cv * inv_diag
return u
end
In [ ]:
function lda(x, y, n_comp)
-- Make y a vector
local y = y:clone():squeeze()
local y_idx = torch.range(1,y:size(1)):long()
local n_classes = y:max()
-- PCA to go to N-C dimensional space
U = pca(x, x:size(1)-n_classes, true)
x = x * U
-- Init
local n = x:size(1)
local d = x:size(2)
-- Mean for each class
local x_mean = torch.zeros(n_classes, d)
local n_class = torch.zeros(n_classes)
for c = 1,n_classes do
local class_idx = y_idx[y:eq(c)]
if class_idx:dim() > 0 then
x_mean[c] = torch.mean(x:index(1, class_idx), 1)
n_class[c] = class_idx:size(1)
end
end
-- Calculate Si and Sw
local Si = torch.Tensor(n_classes, d, d)
local Sw = torch.Tensor(d, d):zero()
for c = 1,n_classes do
local class_idx = y_idx[y:eq(c)]
if class_idx:dim() > 0 then
local x_class = x:index(1, class_idx)
local x_class_m = x_class - torch.ger(torch.ones(x_class:size(1)), x_mean[c]:squeeze())
Si[c] = x_class_m:t() * x_class_m
-- Si[c]:div(n_class[c] - 1)
Sw:add(Si[c])
end
end
-- Calculate S_b
local mean = torch.mean(x, 1)
local Sb = torch.Tensor(d, d):zero()
for c = 1,n_classes do
local mean_class_m = (x_mean[c] - mean):resize(1, d)
local Sb_c = mean_class_m:t() * mean_class_m
Sb_c:mul(n_class[c])
Sb:add(Sb_c)
end
-- Calucalte Sw^{-1}Sb
local Sw_inv = torch.inverse(Sw)
local Sw_inv_Sb = Sw_inv * Sb
-- Get eigenvalues and eigenvectors
local ce, cv = torch.symeig(Sw_inv_Sb, 'V')
-- Sort eigenvalues
local ce, idx = torch.sort(ce, true)
-- Sort eigenvectors
cv = cv:index(2, idx)
-- Keep C-1
ce = ce:sub(1, n_classes-1)
cv = cv:sub(1, -1, 1, n_classes-1)
-- Total transform
W = U * cv
return W
end
In [ ]:
function lpp(x, n_comp)
-- Size
local N = x:size(1)
-- Heat kernel temperature parameter
local t = 1000
-- Epsilon threshold
local eps = 100000
-- Build euclidean distances matrix
local dist = torch.zeros(N, N);
for i = 1,N do
for j = 1,N do
-- Compute norm distance
local dist_i_j = torch.norm(x[{{i}}] - x[{{j}}])
-- If distance < epsilon then not a neighbour
if dist_i_j < eps then
dist[{{i},{j}}] = dist_i_j
else
dist[{{i},{j}}] = 0
end
end
end
-- Build W
local W = torch.exp(-dist / t);
-- Build D, a diagonal matrix continaing the sum of distances of Wi
local D = torch.diag(torch.sum(W, 2):squeeze())
-- Build Laplacian Matrix
local L = D - W
-- LPP
local W_eig = torch.inverse(x:t() * D * x) * x:t() * L * x
-- Solve generalized eigenvector problem to obtain solutions
local ce, cv = torch.symeig(W_eig, 'V')
-- Remove negative and zero eigenvalues
local mask = ce:gt(0)
local mask_idx = torch.range(1,ce:size(1)):long()
ce = ce:index(1, mask_idx[mask])
cv = cv:index(2, mask_idx[mask])
-- Sort eigenvalues in ascending order
local ce, idx = torch.sort(ce)
local cv = cv:index(2, idx)
-- Keep n_comp
if n_comp and n_comp < cv:size(2) then
ce = ce:sub(1, n_comp)
cv = cv:sub(1, -1, 1, n_comp)
end
local A = cv
return A
end
In [ ]:
function npp(x, n_comp)
-- Size
local N = x:size(1)
-- Heat kernel temperature parameter
local t = 1000
-- Epsilon threshold
local eps = 100000
-- Build euclidean distances matrix
local dist = torch.zeros(N, N);
for i = 1,N do
for j = 1,N do
-- Compute norm distance
local dist_i_j = torch.norm(x[{{i}}] - x[{{j}}])
-- If distance < epsilon then not a neighbour
if dist_i_j < eps then
dist[{{i},{j}}] = dist_i_j
else
dist[{{i},{j}}] = 0
end
end
end
-- Build W
local W = torch.exp(-dist / t);
-- Build D, a diagonal matrix continaing the sum of distances of Wi
local D = torch.diag(torch.sum(W, 2):squeeze())
-- Build Laplacian Matrix
local L = D - W
-- NPP
local M = (torch.eye(W:size(1)) - W) * (torch.eye(W:size(1)) - W):t()
local C = x:t() * x
local L = x:t() * M * x
local W_eig = torch.inverse(C) * L
-- Solve generalized eigenvector problem to obtain solutions
local ce, cv = torch.symeig(W_eig, 'V')
-- Remove negative and zero eigenvalues
local mask = ce:gt(0)
local mask_idx = torch.range(1,ce:size(1)):long()
ce = ce:index(1, mask_idx[mask])
cv = cv:index(2, mask_idx[mask])
-- Sort eigenvalues in ascending order
local ce, idx = torch.sort(ce)
local cv = cv:index(2, idx)
-- Keep n_comp
if n_comp and n_comp < cv:size(2) then
ce = ce:sub(1, n_comp)
cv = cv:sub(1, -1, 1, n_comp)
end
local A = cv
return A
end
In [ ]:
function fastica(x, y, n_comp)
-- Functions and Derivatives
local g = torch.tanh
local g_prime = function(u)
return torch.ones(u:size()) - torch.tanh(u):pow(2)
end
-- Init sizes
M = x:size(1)
N = x:size(2)
C = M - 1
-- Center Data
local mean = torch.mean(x, 1)
local x = x - torch.ger(torch.ones(x:size(1)), mean:squeeze())
-- Whiten data
U = pca(x, C, true)
x = x * U
-- Transpose for same notation
x = x:t()
-- Init W and wp with random variables
local W = torch.rand(C, C)
-- Iterate
local max_iter = 100
for p = 1, C do
local converged = false
local iter = 1
-- While Wp changes
while not converged and iter <= max_iter do
local Wp_old = W[{{p}}]:clone()
-- Compute new Wp
local Wp_x = Wp_old * x
local Wp = torch.ones(1, C) * (x * g(Wp_x):t()):mean() - Wp_old * g_prime(Wp_x):mean(2):squeeze()
-- Normalise it
Wp = Wp / torch.norm(Wp)
W[{{p}}] = Wp
local sum_w = torch.zeros(1, C);
if p > 1 then
-- for j = 1, p-1 do
-- sum_w = sum_w + W[{{j}}] * W[{{j}}]:t() * Wp
-- end
sum_w = ((Wp * W[{{1,p-1}}]:t()) * W[{{1,p-1}}])
end
Wp = Wp - sum_w
Wp = Wp / torch.norm(Wp)
-- Check convergence
local lim = 1 - torch.abs(Wp * W[{{p}}]:t()):squeeze()
local converged = lim < 1e-04
iter = iter + 1
if iter > 1 then
W[{{p}}] = Wp
end
end
end
-- De-whitening
ICA = U * W
-- Keep n_comp
if n_comp and n_comp < cv:size(2) then
ICA = ICA:sub(1, -1, 1, n_comp)
end
return ICA
end
In [ ]:
-- Example performing Whitened PCA on X[Samples x Features]
U = pca(X, n_components, true)
X_transformed = X * U
Copyright (C) 2015 John-Alexander Assael (www.johnassael.com) https://github.com/iassael/torch7-decomposition
The MIT License (MIT)
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.