soft.th <- function(lambda, z) {
return(sign(z) * pmax(abs(z) - lambda, 0))
} # z がベクトルでも，soft.th は動作する
SCoTLASS <- function(lambda, X) {
n <- nrow(X)
p <- ncol(X)
v <- rnorm(p)
v <- v / norm(v, "2")
for (k in 1:200) {
u <- X %*% v
u <- u / norm(u, "2")
v <- ## 空欄 ##
v <- soft.th(lambda, v)
size <- norm(v, "2")
if (size > 0) {
v <- v / size
} else {
break
}
}
if (norm(v, "2") == 0)
print("v の全要素が 0 になった")
return(v)
}
## サンプル
n <- 100
p <- 50
X <- matrix(rnorm(n * p), nrow = n)
lambda.seq <- 0:10 / 10
S <- NULL
T <- NULL
for (lambda in lambda.seq) {
v <- SCoTLASS(lambda, X)
S <- c(S, sum(sign(v ^ 2)))
T <- c(T, norm(X %*% v, "2"))
}
plot(lambda.seq, S, xlab = "lambda", ylab = "非ゼロベクトルの個数")
plot(lambda.seq, T, xlab = "lambda", ylab = "分散の和")alphaが不完全なようです。
n <- 100
p <- 10
lambda <- 0.001
mu <- 0.1
x <- matrix(rnorm(n * p), ncol = p)
for (j in 1:p)
x[, j] <- x[, j] - mean(x[, j])
for (j in 1:p)
x[, j] <- x[, j] / sqrt(sum(x[, j] ^ 2))
r <- rep(0, n)
v <- rnorm(p)
for (h in 1:1100) {
z <- x %*% v
theta <- t(x) %*% z / norm(t(x) %*% z, "2")
alpha <- sum(theta ^ 2) +
for (k in 1:p) {
for (i in 1:n)
r[i] <- sum(theta * x[i, ]) - sum(theta ^ 2) * sum(x[i, -k] * v[-k])
S <- ## 空欄(1) ##
T <- ## 空欄(2) ##
v[k] <- soft.th(lambda, S) / T
}
}
v
tk.means <- function(X, K, weights = w) {
n <- nrow(X)
p <- ncol(X)
y <- sample(1:K, n, replace = TRUE)
center <- array(dim = c(K, p))
for (h in 1:10) {
for (k in 1:K) {
if (sum(y[] == k) == 0) {
center[k, ] <- Inf
} else {
for (j in 1:p)
center[k, j] <- ## 空欄(1) ##
}
}
for (i in 1:n) {
S.min <- Inf
for (k in 1:K) {
if (center[k, 1] == Inf)
break
S <- sum((X[i, ] - center[k, ]) ^ 2 * w)
if (S < S.min) {
S.min <- S
## 空欄(2) ##
}
}
}
}
return(y)
}
## データの生成
K <- 10
p <- 2
n <- 1000
X <- matrix(rnorm(p * n), nrow = n, ncol = p)
w <- c(1, 1)
y <- k.means(X, K, w)
## 結果の出力
plot(-3:3, -3:3, xlab = "x", ylab = "y", type = "n")
points(X[, 1], X[, 2], col = y + 1)w.a <- function(a, s) {
if (sum(a) < s)
return(a)
p <- length(a)
lambda <- max(a) / 2
delta <- lambda
for (h in 1:10) {
for (j in 1:p)
w[j] <- soft.th(lambda, ## 空欄(1) ##)
ww <- sqrt(sum(w ^ 2))
if (ww == 0) {
w <- 0
} else {
w <- w / ww
}
if (sum(w) > s) {
lambda <- lambda + delta
} else {
lambda <- ## 空欄(2) ##
}
delta <- delta / 2
}
return(w)
}sparse.k.means <- function(X, K, s) {
p <- ncol(X)
w <- rep(1, p)
for (h in 1:10) {
y <- k.means(## 空欄(1) ##)
a <- comp.a(## 空欄(2) ##)
w <- w.a(## 空欄(3) ##)
}
return(list(w = w, y = y))
}
comp.a <- function(X, y) {
n <- nrow(X)
p <- ncol(X)
a <- array(dim = p)
for (j in 1:p) {
a[j] <- 0
for (i in 1:n) {
for (h in 1:n)
a[j] <- a[j] + (X[i, j] - X[h, j]) ^ 2 / n
}
for (k in 1:K) {
S <- 0
index <- which(y == k)
if (length(index) == 0)
break
for (i in index) {
for (h in index)
S <- S + (X[i, j] - X[h, j]) ^ 2
}
a[j] <- a[j] - S / length(index)
}
}
return(a)
}
## 実行
p <- 10
n <- 100
X <- matrix(rnorm(p * n), nrow = n, ncol = p)
sparse.k.means(X, 5, 1.5)ww <- function(x, mu = 1, dd = 0) {
n <- nrow(x)
w <- array(dim = c(n, n))
for (i in 1:n) {
for (j in 1:n)
w[i, j] <- exp(-mu * sum((x[i, ] - x[j, ]) ^ 2))
}
if (dd > 0) {
for (i in 1:n) {
dis <- NULL
for (j in 1:n)
dis <- c(dis, sqrt(sum((x[i, ] - x[j, ]) ^ 2)))
index <- which(dis > dd)
w[i, index] <- 0
}
}
return(w)
}
prox <- function(x, tau) {
if (sum(x ^ 2) == 0) {
return(x)
} else {
return(max(0, 1 - tau / sqrt(sum(x ^ 2))) * x)
}
}
update.u <- function(v, lambda) {
u <- array(dim = c(n, d))
z <- 0
for (i in 1:n)
z <- z + x[i, ]
y <- x
for (i in 1:n) {
if (i < n) {
for (j in (i + 1):n)
y[i, ] <- y[i, ] + lambda[i, j, ] + nu * v[i, j, ]
}
if (1 < i) {
for (j in 1:(i - 1))
y[i, ] <- y[i, ] - lambda[j, i, ] - nu * v[j, i, ]
}
u[i, ] <- (y[i, ] + nu * z) / (n * nu + 1)
}
return(u)
}
update.v <- function(u, lambda) {
v <- array(dim = c(n, n, d))
for (i in 1:(n - 1)) {
for (j in (i + 1):n)
v[i, j, ] <- prox(u[i, ] - u[j, ] - lambda[i, j, ] / nu, gamma * w[i, j] / nu)
}
return(v)
}
update.lambda <- function(u, v, lambda) {
for (i in 1:(n - 1)) {
for (j in (i + 1):n)
lambda[i, j, ] <- lambda[i, j, ] + nu * (v[i, j, ] - u[i, ] + u[j, ])
}
return(lambda)
}
## u, v, lambda の更新のサイクルを max_iter だけ繰り返す
convex.cluster <- function() {
v <- array(rnorm(n * n * d), dim = c(n, n, d))
lambda <- array(rnorm(n * n * d), dim = c(n, n, d))
for (iter in 1:max_iter) {
u <- ## 空欄(1) ##
v <- ## 空欄(2) ##
lambda <- ## 空欄(3) ##
}
return(list(u = u, v = v))
}
## データ生成
n <- 50
d <- 2
x <- matrix(rnorm(n * d), n, d)
## 凸クラスタリングの実行
w <- ww(x, 1, dd = 1.5)
gamma <- 50
nu <- 1
max_iter <- 100
v <- convex.cluster()$v
## 隣接行列の計算
a <- array(0, dim = c(n, n))
for (i in 1:(n - 1)) {
for (j in (i + 1):n)
if (sqrt(sum(v[i, j, ] ^ 2)) < 1 / 10 ^ 4) {
a[i, j] <- 1
a[j, i] <- 1
}
}
}
## 作図
k <- 0
y <- rep(0, n)
for (i in 1:n) {
if (y[i] == 0) {
k <- k + 1
y[i] <- k
if (i < n) {
for (j in (i + 1):n) {
if (a[i, j] == 1)
y[j] <- k
}
}
}
}
plot(0, xlim = c(-3, 3), ylim = c(-3, 3), type = "n", main = "gamma = 50")
points(x[, 1], x[, 2], col = y + 1)s.update.u <- function(v, lambda) {
u <- array(dim = c(n, d))
y <- x
for (i in 1:n) {
if (i < n) {
for (j in (i + 1):n)
y[i, ] <- y[i, ] + lambda[i, j, ] + nu * v[i, j, ]
}
if (1 < i) {
for (j in 1:(i - 1))
y[i, ] <- y[i, ] - lambda[j, i, ] - nu * v[j, i, ]
}
}
for (j in 1:d)
u[, j] <- gr(G, G.inv %*% y[, j], gamma.2 * r[j])
for (j in 1:d)
u[, j] <- u[, j] - mean(u[, j])
return(u)
}
s.convex.cluster <- function() {
## gamma.2, r[1], ..., r[p] を設定しておく
G <- sqrt(1 + n * nu) * diag(n) - (sqrt(1 + n * nu) - 1) / n %*% matrix(1, n, n)
G.inv <- (1 + n * nu) ^ (-0.5) %*% (diag(n) + (sqrt(1 + n * nu) - 1) / n *
matrix(1, n, n))
v <- array(rnorm(n * n * d), dim = c(n, n, d))
lambda <- array(rnorm(n * n * d), dim = c(n, n, d))
for (iter in 1:max_iter) {
u <- s.update.u(v, lambda)
v <- update.v(u, lambda)
lambda <- update.lambda(u, v, lambda)
}
return(list(u = u, v = v))
}