[
https://issues.apache.org/jira/browse/SYSTEMML-510?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel
]
Mike Dusenberry updated SYSTEMML-510:
-------------------------------------
Description:
If we look at the inner loop of Poisson nonnegative matrix factorization (PNMF)
in general, we update the factors as
{code}
H = (H * (t(W) %*% (V/(W%*%H + 1e-17))))/t(colSums(W)
W = (W * ((V/(W%*%H + 1e-17)) %*% t(H)))/t(rowSums(H))
{code}.
Notice the addition of the "1e-17" epsilon term in the denominators.
Mathematically, we need this in case any cell of {{W%*%H}} evaluates to zero so
that we can avoid dividing by zero. R needs this, but SystemML technically
does not due to a fused operator, "wdivmm", that takes care of these situations
(or this may be done in the general case?). This fused operator is currently
applied to the pattern {{t(W) %*% (V / %* (W %*% H))}}, amongst other similar
patterns. Ideally, this would easily apply to {{t(W) %*% (V/(W%*%H + 1e-17)}},
regardless of the unneeded epsilon term. Currently, the addition of the
epsilon term causes the algorithm to run in non-linear time (quad or
exponential). Initially, the behavior pointed towards the possibility of the
optimizer avoiding the rewrite to the fused operator, resulting in naive
computation, and non-linear growth in training time. Further exploration seems
to show that the rewrite is indeed still being applied, but there seems to also
be a recursive nesting of the same rewrite over various regions of the above
statements that is not found when the epsilon term is removed.
The following is the full PNMF DML script used:
{code}
V = read($X)
max_iteration = $maxiter
rank = $rank
n = nrow(V)
m = ncol(V)
range = 0.01
W = Rand(rows=n, cols=rank, min=0, max=range, pdf="uniform")
H = Rand(rows=rank, cols=m, min=0, max=range, pdf="uniform")
loglik0 = sum(V*log(W%*%H)) - as.scalar(colSums(W)%*%rowSums(H))
i=0
while(i < max_iteration) {
# Addition of epsilon (1e-17) term causes script to run in non-linear time:
H = (H * (t(W) %*% (V/(W%*%H + 1e-17))))/t(colSums(W))
W = (W * ((V/(W%*%H + 1e-17)) %*% t(H)))/t(rowSums(H))
# Removal of epsilon works correctly:
#H = (H * (t(W) %*% (V/(W%*%H))))/t(colSums(W))
#W = (W * ((V/(W%*%H)) %*% t(H)))/t(rowSums(H))
i = i + 1;
print(i + "")
}
loglik = sum(V*log(W%*%H+1e-17)) - as.scalar(colSums(W)%*%rowSums(H))
print("pnmf: " + loglik0 + " -> " + loglik)
{code}
was:
If we look at the inner loop of Poisson nonnegative matrix factorization (PNMF)
in general, we update the factors as
{code}
H = (H * (t(W) %*% (V/(W%*%H + 1e-17))))/t(colSums(W)
W = (W * ((V/(W%*%H + 1e-17)) %*% t(H)))/t(rowSums(H))
{code}.
Notice the addition of the "1e-17" epsilon term in the denominators.
Mathematically, we need this in case any cell of W%*%H evaluates to zero so
that we can avoid dividing by zero. R needs this, but SystemML technically
does not due to a fused operator, "wdivmm", that takes care of these situations
(or this may be done in the general case?). This fused operator is currently
applied to the pattern {{t(W) %*% (V / %* (W %*% H))}}, amongst other similar
patterns. Ideally, this would easily apply to {{t(W) %*% (V/(W%*%H + 1e-17)}},
regardless of the unneeded epsilon term. Currently, the addition of the
epsilon term causes the algorithm to run in non-linear time (quad or
exponential). Initially, the behavior pointed towards the possibility of the
optimizer avoiding the rewrite to the fused operator, resulting in naive
computation, and non-linear growth in training time. Further exploration seems
to show that the rewrite is indeed still being applied, but there seems to also
be a recursive nesting of the same rewrite over various regions of the above
statements that is not found when the epsilon term is removed.
The following is the full PNMF DML script used:
{code}
V = read($X)
max_iteration = $maxiter
rank = $rank
n = nrow(V)
m = ncol(V)
range = 0.01
W = Rand(rows=n, cols=rank, min=0, max=range, pdf="uniform")
H = Rand(rows=rank, cols=m, min=0, max=range, pdf="uniform")
loglik0 = sum(V*log(W%*%H)) - as.scalar(colSums(W)%*%rowSums(H))
i=0
while(i < max_iteration) {
# Addition of epsilon (1e-17) term causes script to run in non-linear time:
H = (H * (t(W) %*% (V/(W%*%H + 1e-17))))/t(colSums(W))
W = (W * ((V/(W%*%H + 1e-17)) %*% t(H)))/t(rowSums(H))
# Removal of epsilon works correctly:
#H = (H * (t(W) %*% (V/(W%*%H))))/t(colSums(W))
#W = (W * ((V/(W%*%H)) %*% t(H)))/t(rowSums(H))
i = i + 1;
print(i + "")
}
loglik = sum(V*log(W%*%H+1e-17)) - as.scalar(colSums(W)%*%rowSums(H))
print("pnmf: " + loglik0 + " -> " + loglik)
{code}
> Generalize All Existing wdivmm Patterns
> ---------------------------------------
>
> Key: SYSTEMML-510
> URL: https://issues.apache.org/jira/browse/SYSTEMML-510
> Project: SystemML
> Issue Type: Bug
> Components: Compiler, Parser, Runtime
> Reporter: Mike Dusenberry
>
> If we look at the inner loop of Poisson nonnegative matrix factorization
> (PNMF) in general, we update the factors as
> {code}
> H = (H * (t(W) %*% (V/(W%*%H + 1e-17))))/t(colSums(W)
> W = (W * ((V/(W%*%H + 1e-17)) %*% t(H)))/t(rowSums(H))
> {code}.
> Notice the addition of the "1e-17" epsilon term in the denominators.
> Mathematically, we need this in case any cell of {{W%*%H}} evaluates to zero
> so that we can avoid dividing by zero. R needs this, but SystemML
> technically does not due to a fused operator, "wdivmm", that takes care of
> these situations (or this may be done in the general case?). This fused
> operator is currently applied to the pattern {{t(W) %*% (V / %* (W %*% H))}},
> amongst other similar patterns. Ideally, this would easily apply to {{t(W)
> %*% (V/(W%*%H + 1e-17)}}, regardless of the unneeded epsilon term.
> Currently, the addition of the epsilon term causes the algorithm to run in
> non-linear time (quad or exponential). Initially, the behavior pointed
> towards the possibility of the optimizer avoiding the rewrite to the fused
> operator, resulting in naive computation, and non-linear growth in training
> time. Further exploration seems to show that the rewrite is indeed still
> being applied, but there seems to also be a recursive nesting of the same
> rewrite over various regions of the above statements that is not found when
> the epsilon term is removed.
> The following is the full PNMF DML script used:
> {code}
> V = read($X)
> max_iteration = $maxiter
> rank = $rank
> n = nrow(V)
> m = ncol(V)
> range = 0.01
> W = Rand(rows=n, cols=rank, min=0, max=range, pdf="uniform")
> H = Rand(rows=rank, cols=m, min=0, max=range, pdf="uniform")
> loglik0 = sum(V*log(W%*%H)) - as.scalar(colSums(W)%*%rowSums(H))
> i=0
> while(i < max_iteration) {
> # Addition of epsilon (1e-17) term causes script to run in non-linear time:
> H = (H * (t(W) %*% (V/(W%*%H + 1e-17))))/t(colSums(W))
> W = (W * ((V/(W%*%H + 1e-17)) %*% t(H)))/t(rowSums(H))
> # Removal of epsilon works correctly:
> #H = (H * (t(W) %*% (V/(W%*%H))))/t(colSums(W))
> #W = (W * ((V/(W%*%H)) %*% t(H)))/t(rowSums(H))
> i = i + 1;
> print(i + "")
> }
> loglik = sum(V*log(W%*%H+1e-17)) - as.scalar(colSums(W)%*%rowSums(H))
> print("pnmf: " + loglik0 + " -> " + loglik)
> {code}
--
This message was sent by Atlassian JIRA
(v6.3.4#6332)
