Assess the dependence of LDA on the random seed New to latent Dirichlet Allocation (LDA), I would like to be sure that my output (in the first the step, the word-per-topic probabilities) depends on the input merely, and is (somewhat) stable whatever start value for the first assignment is used. 
How can I assess that? I assume this is not in-built in the algorithm (I am using the R packages topicmodels and lda).
More specifically, assume the number of topics is 2 and I run the LDA 3 times with identical input but with different random seeds. Now the word-per-topic distribution might look like:
# run 1
list(
  topic1=c(word1=0.1, word2=0.05, word3=0.01, ...),
  topic2=c(word4=0.2, word5=0.15, word6=0.1, ...)
)
# run 2
list(
  topic1=c(word1=0.2, word7=0.05, word3=0.01, ...),
  topic2=c(word4=0.2, word5=0.15, word2=0.05, ...)
)
# run 3
list(
  topic1=c(word4=0.2, word5=0.15, word9=0.05, ...),
  topic2=c(word3=0.1, word1=0.05, word2=0.01, ...)
)

Each element of the list (topic1, topic2) is a numeric vector on the whole vocabulary that sums to 1.
How "similar" are the outputs of these three runs? In particular, I don't want to penalize different orders of the results (it is IMHO natural that sometimes the words are assigned to the first topic, and, with a different start seed, the same words to the second topic).
EDIT: As an example, consider the following LDA (using R):
library(lda)
data(cora.documents)
data(cora.vocab)
K <- 10 ## Num clusters

one_run <- function(seed) {
    set.seed(seed)
    fit <- lda.collapsed.gibbs.sampler(cora.documents, 
        K,  ## Num clusters
        cora.vocab,
        25,  ## Num iterations
        0.1,
        0.1,
        compute.log.likelihood=FALSE
    ) 
    return(fit$topics/(rowSums(fit$topics) + 1e-05))   
}

Here, the function one_run performs the LDA with fixed input except for the random start value. Now I run this function 3 times:
seeds <- c(123, 456, 789)
res <- lapply(seeds, one_run)

and I get a list of the word-topic distribution for each run. For example,
res[[1]] is a matrix with 10 rows (topics) and 2961 columns (for each term of the vocabulary one column) that contains the probability that term j belongs to topic i. Now I would like to assess the dependence of the output given the random seed; however the order of the rows might differ. For instance, if I look at the words with most of the probability mass:
lapply(res, top.topic.words, num.words=5)

it appears that topic 1 from the first run (with the top 5 words "learning", "problem", "paper", "system", "control") corresponds to topic 2 from the second run (with the words "learning", "problem", "system", "paper", "theory"). I want to assess the dependence from the random seed, but ignore the order of the topics in the different runs.
What is the best practice here? I first was thinking a modified (probability-weighted) version of the Jaccard-Index, but I don't know how to deal with the order issue. Now I am thinking maybe some kind of assignment problem, but I am still hoping for an easier (faster) way.
I think the problem is similar when assessing the output of different runs of the k-means algorithm.
 A: Here is what I came up with. Maybe it is still buggy, certainly there is room for improvement.
I programmed a linear assignment problem: given the outputs of two different runs and a distance function, it reorders one of the outputs so that the difference is minimal. This minimal difference is returned.
library(lpSolveAPI)
similarity_order <- function(a, b, distance, browse=FALSE) {
    if(!is.matrix(a) && !is.matrix(b)) {
        a <- matrix(a, ncol=1)
        b <- matrix(b, ncol=1)
    }
    n <- nrow(a)

    # build distance matrix
    dist_matrix <- matrix(0, nrow=n, ncol=n)
    for (i in seq(1, n))
        for (j in seq(1, n)) {
            dist_matrix[i,j] <- distance(a[i,], b[j,])
            # if(j>i) dist_matrix[j,i] <- dist_matrix[i,j] # use symmetry
        }

    # setup and solve linear problem
    lprec <- make.lp(nrow=2*n, ncol=n^2) # z,s -> (s-1)*n+z
    set.objfn(lprec, as.vector(dist_matrix))
    set.type(lprec, seq(1,n^2), "binary")
    for (k in seq(1,n)) {
        # column sums == 1
        set.row(lprec, k, xt=rep(1,n), indices=seq((k-1)*n+1, length.out=n, by=1))
        # row sums == 1
        set.row(lprec, n+k, xt=rep(1,n), indices=k+seq(0, length.out=n, by=n))
    }
    set.rhs(lprec, b=rep(1,2*n))
    set.constr.type(lprec, types=rep("=", 2*n))
    solve(lprec)

    # calculate the distance between the reordered a and b and return its sum
    ans <- get.variables(lprec)
    idx <- apply(matrix(ans, nrow=n),2,function(x) which(x==1))
    ret <- rep(NA, n)
    for (k in seq(1,n)) ret[k] <- distance(a[idx[k],], b[k,]) # maybe that is just get.objective(lprec)?
    if(browse) browser()
    return(sum(ret))
}

To assess the dependence on the seeds, the following wrapper function just calls the algorithm repeatedly and calculates the minimal distances over all combinations of the outputs.
seed_stability <- function(fun, seeds, distance, ...) {
    if(length(seeds)==1) seeds <- sample.int(n=.Machine$integer.max, seeds, replace=FALSE)
    fits <- lapply(seeds, function(x) {set.seed(x); fun(...)})
    arg <- combn(seq(1, length(seeds)), 2)
    distances <- rep(NA, ncol(arg))
    for (i in seq(1, ncol(arg))) distances[i] <- similarity_order(fits[[arg[1,i]]], fits[[arg[2,i]]], distance=distance)
    # browser()
    return(sum(distances)/length(distances))
}

Now to assess a LDA procedure, I use the Hellinger distance function:
hellinger_distance <- function(x,y) sqrt(sum((sqrt(x)-sqrt(y))^2)/2)

For the example in the question:
test_cora <- function(num.iterations=50) {
    data(cora.documents, package="lda")
    data(cora.vocab, package="lda")
    K <- 10 ## Num clusters

    one_run <- function() {
        fit <- lda.collapsed.gibbs.sampler(cora.documents, 
            K,  ## Num clusters
            cora.vocab,
            num.iterations=num.iterations,  ## Num iterations
            alpha=0.1,
            eta=0.1,
            compute.log.likelihood=FALSE
        ) 
        return(fit$topics/(rowSums(fit$topics) + 1e-05))   
    }
    seed_stability(one_run, seeds=20, distance=hellinger_distance)
}

For example, varying the number of iterations parameter, the test_cora shows that the dependence on the seeds reduces from 5.9 to 5.3 when the number of iterations is increased from 25 to 50.
The function seed_stability can also be used to assess the output of kmeans. Here I use as distance function the absolute norm, 
abs_distance <- function(x,y) sum(abs(x-y))

and then the following function tests kmeans on two datasets:
test_kmeans <- function(good=TRUE) {
    x <- c(0,10)
    y <- c(20,15)
    z <- c(10,8)
    nice_dat <- cbind(
        c(x[1]+rnorm(5), y[1]+rnorm(5), z[1]+rnorm(5)),
        c(x[2]+rnorm(5), y[2]+rnorm(5), z[2]+rnorm(5))
    )

    bad_dat <- cbind(runif(50),runif(50))       

    one_run <- function(dat, k) {
        fit <- kmeans(dat, centers=k, nstart=3)
        plot(x=dat[,1], y=dat[,2], col=fit$cluster)
   return(ifelse(t((sapply(1:k, function(x) fit$cluster==x))), 1, 0))
    }

    par(mfrow=c(2,3))
    if(good) 
        seed_stability(one_run, seeds=6, dat=nice_dat, k=3, distance=abs_distance)
    else
        seed_stability(one_run, seeds=6, dat=bad_dat, k=3, distance=abs_distance)
}

It shows that for the nicely clustered dataset the algorithm returns always 0, hence there is no dependence on the seeds. For the uniform distributed dataset, however, there is often a dependency found.
For now, I am happy with that approach, but I am still looking for other ideas.
