Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commit6ccddfd

Browse files
committed
add tucker.nonneg()
add nonnegative Tucker decomposition
1 parenta9d607e commit6ccddfd

File tree

2 files changed

+207
-0
lines changed

2 files changed

+207
-0
lines changed

‎rTensor/NAMESPACE‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ export(tperm)
2929
export(ttl)
3030
export(ttm)
3131
export(tucker)
32+
export(tucker.nonneg)
3233
export(unfold)
3334
export(unmatvec)
3435
export(vec)

‎rTensor/R/rTensor_Decomp.R‎

Lines changed: 206 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -516,3 +516,209 @@ t_svd_reconstruct <- function(L){
516516
resid<- fnorm(tnsr-est)
517517
invisible(list(core= as.tensor(core),est=est,fnorm_resid=resid,norm_percent= (1-resid/fnorm(tnsr))*100))
518518
}
519+
520+
#' nonnegative Tucker decomposition
521+
#'
522+
#' Decomposes nonnegative tensor \code{tnsr} into core nonnegative tensor \code{Z} and nonnegative factor matrices \code{U[n]}.
523+
#'@export
524+
#'@param tnsr nonnegative tensor with \code{K} modes
525+
#'@param ranks an integer vector of length \code{K} specifying the modes sizes for the output core tensor \code{Z}
526+
#'@param tol relative Frobenius norm error tolerance
527+
#'@param max_iter maximum number of iterations if error stays above \code{tol}
528+
#'@param max_time max running time
529+
#'@param rw controls the extrapolation weight
530+
#'@param U0 initial factor matrices, defaults to nonnegative Gaussian random matrices
531+
#'@param Z0 initial core tensor \code{Z}, defaults to nonnegative Gaussian random tensor
532+
#'@param verbose more output algorithm progress
533+
#'@return a list:\describe{
534+
#'\item{\code{U}}{nonnegative factor matrices}
535+
#'\item{\code{Z}}{nonnegative core tensor}
536+
#'\item{\code{est}}{estimate \eqn{Z \times_1 U_1 \ldots \times_K U_K}}
537+
#'\item{\code{conv}}{method convergence indicator}
538+
#'\item{\code{fnorm_resid}}{the Frobenius norm of the error \code{fnorm(est-tnsr)} - if there was no truncation, then this is O(mach_eps) }
539+
#'\item{\code{n_iter}}{number of iterations}
540+
#'\item{\code{diag}}{convergence info for each iteration\describe{
541+
#'\item{\code{all_resids}}{residues}
542+
#'\item{\code{all_rel_resid_deltas}}{residue delta relative to the current residue}
543+
#'\item{\code{all_rel_resids}}{residue relative to the \code{sqrt(||tnsr||)}}
544+
#'}}}
545+
#'
546+
#'@details The algorithm uses the block-coordinate update procedure to solve the following optimization problem:
547+
#' \deqn{\min 0.5 \|tnsr - Z \times_1 U_1 \ldots \times_K U_K \|_{F^2}, \;\text{where}\; Z \geq 0, \, U_i \geq 0.}
548+
#' It stops if either the relative improvement of the error is below the tolerance \code{tol} for 3 consequitive iterations or
549+
#' both the relative error improvement and relative error (wrt the \code{tnsr} norm) are below the tolerance.
550+
#' Otherwise it stops if the maximal number of iterations or the time limit were reached.
551+
#'
552+
#'@note The implementation is based on ntd() MATLAB code by Yangyang Xu and Wotao Yin.
553+
#'@references Y. Xu and W. Yin, "A Block Coordinate Descent Method for Multi-Convex Optimization with Applications to Nonnegative Tensor Factorization and Completion". SIAM Journal on Imaging Sciences, 6(3), 1758-1789, 2013.
554+
#'@seealso \code{\link{tucker}}
555+
#'@seealso \url{http://www.math.ucla.edu/~wotaoyin/papers/bcu/ntd/ntd.html}
556+
#'@seealso \url{http://www.math.ucla.edu/~wotaoyin/papers/bcu/index.html}
557+
tucker.nonneg<-function(tnsr,ranks,
558+
tol=1e-4,
559+
max_iter=500,max_time=0,rw=1,
560+
U0=NULL,Z0=NULL,verbose=FALSE )
561+
{
562+
make_nonneg.tnsr<-function(tnsr )
563+
{
564+
tnsr@data[tnsr@data<0 ]<-0
565+
return (tnsr )
566+
}
567+
make_nonneg.mtx<-function(mtx )
568+
{
569+
mtx[mtx<0 ]<-0
570+
return (mtx )
571+
}
572+
Kway<- dim(tnsr)# dimension of tnsr
573+
K<- length(Kway)# tnsr is an K-way tensor
574+
575+
if ( is.null(U0) ) {
576+
if (verbose ) message('Generating random initial factor matrices estimates...' )
577+
U0<- lapply( seq_len(K),function(n) make_nonneg.mtx(matrix( rnorm(Kway[[n]]*ranks[[n]] ),ncol=ranks[[n]] ) ) )
578+
}
579+
if ( is.null(Z0) ) {
580+
if (verbose ) message('Generating random initial core tensor estimate...' )
581+
Z0<- make_nonneg.tnsr( rand_tensor(modes=ranks,drop=FALSE) )
582+
}
583+
584+
Tnrm<- fnorm(tnsr)
585+
fnorm_resid0<-0.5*Tnrm^2
586+
587+
U0<- lapply(U0,function(U0i)U0i/norm(U0i,"F")*Tnrm^(1/(K+1)) )
588+
Usq<- lapply(U0,crossprod )
589+
Z0<-Z0/fnorm(Z0)*Tnrm^(1/(K+1))
590+
591+
U<-U0
592+
Um<-U0
593+
Z<-Z0
594+
Zm<-Z0
595+
596+
n_stall<-0
597+
w<-0
598+
t0<-1
599+
t<-t0
600+
wU<- rep.int(1,K+1 )
601+
L0<-wU
602+
L<-L0
603+
604+
# Precompute matrix unfoldings of input tensor to save computing time if it is not too large
605+
if (K*prod(Kway)<4000^2) {
606+
if (verbose ) message('Precomputing input tensor unfoldings...' )
607+
Tmtx<- lapply( seq_len(K),function(n) unfold(tnsr,row_idx=n,col_idx= seq_len(K)[-n] )@data )
608+
}else {
609+
if (verbose ) message('Input tensor too big, no unfoldings precomputation...' )
610+
Tmtx<-NULL
611+
}
612+
613+
# Iterations of block-coordinate update
614+
# iteratively updated variables:
615+
# GradU: gradients with respect to each component matrix of U
616+
# GradZ: gradient with respect to Z
617+
# U,Z: new updates
618+
# U0,Z0: old updates
619+
# Um,Zm: extrapolations of U
620+
# L, L0: current and previous Lipschitz bounds
621+
# resid, resid0: current and previous residual error
622+
start_time<- proc.time()
623+
#progress bar
624+
pb<- txtProgressBar(min=0,max=max_iter,style=3)
625+
all_resids<-numeric(0)
626+
all_rel_resid_deltas<-numeric(0)
627+
all_rel_resids<-numeric(0)
628+
629+
conv<-FALSE
630+
for (n_iterin seq_len(max_iter)) {
631+
setTxtProgressBar(pb,n_iter)
632+
633+
# -- update the core tensor Z --
634+
L0[[K+1]]<-L[[K+1]]
635+
L[[K+1]]<- prod( sapply(Usq,norm,type='F' ) )
636+
637+
# compute the gradient
638+
GradZ<- ttl(Zm,Usq, seq_len(K))- ttl(tnsr,U, seq_len(K),transpose=TRUE )
639+
# update core tensor
640+
Z<- make_nonneg.tnsr(Zm-GradZ/L[[K+1]] )
641+
642+
# -- update factor matrices U --
643+
for (nin seq_len(K)) {
644+
if (!is.null(Tmtx)) {
645+
B<- unfold( ttl(Z,U[-n], seq_len(K)[-n] ),n, seq_len(K)[-n] )
646+
Bsq<- tcrossprod(B@data,B@data)
647+
MB<- tcrossprod(Tmtx[[n]],B@data)
648+
}else {
649+
Z.mtx<- unfold(Z,n, seq_len(K)[-n])
650+
Bsq<- unfold( ttl(Z,Usq[-n], seq_len(K)[-n] ),n, seq_len(K)[-n] )
651+
Bsq<- tcrossprod(Bsq@data,Z.mtx@data )
652+
653+
MB<- unfold( ttl(tnsr,U[-n], seq_len(K)[-n],transpose=TRUE ),n, seq_len(K)[-n] )
654+
MB<- tcrossprod(MB@data,Z.mtx@data )
655+
}
656+
#compute the gradient
657+
GradU<-Um[[n]]%*%Bsq-MB
658+
L0[[n]]<-L[[n]]
659+
L[[n]]<- norm(Bsq,'F')
660+
U[[n]]<- make_nonneg.mtx(Um[[n]]-GradU/L[[n]] )
661+
Usq[[n]]<- crossprod(U[[n]],U[[n]] )
662+
}
663+
664+
# --- diagnostics, reporting, stopping checks ---
665+
fnorm_resid<- max(0,0.5*(sum(Usq[[K]]*Bsq)-2*sum(U[[K]]*MB)+Tnrm^2) )
666+
rel_resid_delta<- abs(fnorm_resid-fnorm_resid0)/(fnorm_resid0+1)
667+
rel_resid<- sqrt(2*fnorm_resid)/Tnrm
668+
669+
# reporting
670+
all_resids<- append(all_resids,fnorm_resid )
671+
all_rel_resid_deltas<- append(all_rel_resid_deltas,rel_resid_delta )
672+
all_rel_resids<- append(all_rel_resids,rel_resid )
673+
674+
# check stopping criterion
675+
crit<-rel_resid_delta<tol
676+
n_stall<- ifelse(crit,n_stall+1,0 )
677+
if (n_stall>=3||rel_resid<tol ) {
678+
if (verbose ) {
679+
if (rel_resid==0 ) message('Residue is zero. Exact decomposition was found' )
680+
if (n_stall>=3 ) message('Residue relative delta below',tol,'',n_stall,' times in a row' )
681+
if (rel_resid<tol ) message('Residue is',rel_resid,' times below input tensor norm' )
682+
message('tucker.nonneg() converged in',n_iter,' iteration(s),',n_redo,' redo steps' )
683+
}
684+
conv=TRUE
685+
break
686+
}
687+
if (max_time>0&& ( proc.time()-start_time )[[3]]>max_time ) {
688+
warning("Maximal time exceeded, might be not an optimal solution")
689+
break
690+
}
691+
692+
# --- correction and extrapolation ---
693+
t<- (1+sqrt(1+4*t0^2))/2
694+
if (fnorm_resid>=fnorm_resid0) {
695+
# restore U, Z to make the objective nonincreasing
696+
Um<-U0
697+
Zm<-Z0
698+
}else {
699+
#get new extrapolated points
700+
w<- (t0-1)/t# extrapolation weight
701+
# choose smaller weight for convergence
702+
wU<- pmin(w,rw*sqrt(L0/L) )
703+
Um<- lapply( seq_len(K),function(n)U[[n]]+wU[[n]]*(U[[n]]-U0[[n]]) )
704+
Zm<-Z+wU[K+1]*(Z-Z0)
705+
U0<-U
706+
Z0<-Z
707+
t0<-t
708+
fnorm_resid0<-fnorm_resid
709+
}
710+
}
711+
setTxtProgressBar(pb,max_iter)
712+
close(pb)
713+
if (!conv&&n_iter==max_iter ) {
714+
warning("Maximal number of iterations reached, might be not an optimal solution")
715+
}
716+
717+
return (invisible(list(U=U,Z=Z,est=ttl(Z,U,seq_len(K)),
718+
n_iter=n_iter,
719+
conv=conv,
720+
fnorm_resid=fnorm_resid,
721+
diag=list(all_resids=all_resids,
722+
all_rel_resid_deltas=all_rel_resid_deltas,
723+
all_rel_resids=all_rel_resids ) ) ) )
724+
}

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp