Как оптимизировать поиск следа умножения квадратной матрицы?

Question

Я пытаюсь оптимизировать функцию spdep для R для моего случая использования, поскольку она очень медленная для больших баз данных. У меня все было в основном хорошо, но я застрял в одном месте, где я пытаюсь найти след моей матрицы весов для теста ошибок LM. Я думаю, что формула tr[(W' + W) W] (стр. 82 Анселин Л., Бера А.К., Флоракс Р. и Юн М.Дж., 1996 г. Простые диагностические тесты пространственной зависимости. Региональная наука и городская экономика, 26, 77–104.) W - квадратная матрица весов, содержащая пространственную связь каждого наблюдения с другим. Операция tr() - это сумма диагоналей.

В моем случае матрица весов симметрична, а диагонали равны нулю. Итак, я думал, что формула tr [(W '+ W) W] равна 2*sumsq(W), что очень быстро. Но, очевидно, я где-то ошибаюсь, потому что результаты не совпадают с результатами библиотеки spdep, что, вероятно, будет правильным.

Соответствующая часть библиотеки spdep находится здесь. Кто-нибудь может мне помочь, чем результат следующей функции отличается от 2 * sumsq (W) или как сделать это намного быстрее? В этой функции функция lm.LMtests забивается для больших наборов данных.

tracew <- function (listw) {
    dlmtr <- 0
    n <- length(listw$neighbours)
    if (n < 1) stop("non-positive n")
    ndij <- card(listw$neighbours)
    dlmtr <- 0
    for (i in 1:n) {
        dij <- listw$neighbours[[i]]
        wdij <- listw$weights[[i]]
        for (j in seq(length=ndij[i])) {
            k <- dij[j]
# Luc Anselin 2006-11-11 problem with asymmetric listw
                dk <- which(listw$neighbours[[k]] == i)
                if (length(dk) > 0L && dk > 0L &&
                dk <= length(listw$neighbours[[k]]))
                wdk <- listw$weights[[k]][dk]
                else wdk <- 0
                dlmtr <- dlmtr + (wdij[j]*wdij[j]) + (wdij[j]*wdk)
        }
    }
    dlmtr
}

Дополнительное объяснение для тех, кто не знаком с библиотекой spdep языка R: на входе функции listw содержится "графическая" реализация весовой матрицы с двумя списками списков. listw$ соседей - это список, где каждый элемент списка представляет собой список индексов наблюдений, к которым относится наблюдение. listw$ взвешивает список той же структуры с соседями, за исключением того, что он содержит веса отношения.

Заранее спасибо за любые комментарии и указания.

# example code

# initiliaze
library(spdep)
library(multiway)
# load the tracew function above
data(columbus)
columbus = columbus[rep(row.names(columbus), 20), ] # the difference becomes dramatic when n is high. try not replicating at first to see the results.

# manual calculation, using sumsq
w = distm(cbind(columbus$X, columbus$Y))
w[w > 1000000] = Inf # remove some relations acc. to pre-defined rule
w = 1/(1+w)
diag(w) = 0
w = w / (sum(w) / length(columbus$X)) #"C style" standardization
2*sumsq(w)

# spdep calculation
neighs.band = dnearneigh(cbind(columbus$X, columbus$Y), 0, 1000, longlat = TRUE)
w.spdep = lapply(nbdists(neighs.band, cbind(columbus$X, columbus$Y), longlat = TRUE), function(x) 1/(0.001+x))
my.listw = nb2listw(neighs.band, glist = w.spdep, style="C")
tracew(my.listw)