Применение гравитационной модели для оценки матрицы отключения в r

Question

Я применяю гравитационную модель для оценки матрицы двойных ограничений, в R.

Эта модель говорит, что число поездок между отправлением i и назначением j определяется как:

T_ {i} = A_ {i} O_ {i} B_ {j} D_ {j} f (c_ {ij})

Где:

A_ {i} = \ left [\ sum_ {j} {B_ {j} D_ {j} f (c_ {ij})} \ right] ^ {- 1}

и

B_ {j} = \ left[\ sum_ {i} {A_ {i} O_ {i} f (c_ {ij})} \ right] ^ {- 1}

A_ {i} и B_ {j} являются параметрами, которые должныоценивать так, чтобы: \ sum_ {j} {T_ {ij}} = O_ {i} и \ sum_ {i} {T_ {ij}} = D_ {i}

Кроме того, в этом случае Iустановит:

f (c_ {ij}) = c ^ {- 0,9}

Это просто мера стоимости поездки («импеданс») между парами источник-пункт назначения.

Матрица отключения по этой модели может быть оценена с помощью процедуры итеративного двупропорционального подбора («поставки»), которая неявно вычисляет A_ {i} и B_ {j}.

Следующий код делает это для воображаемого набора итогов отправления и назначения и некоторых случайно выбранных данных о затратах (где стоимость должна быть строго положительной) для матрицы 5 на 5.

Основная идея состоит в том, чтобы сначала ограничить суммы строк матрицы, чтобы они равнялись заданным итоговым значениям исходной поездки, а затем ограничить суммы столбцов матрицы равными заданным суммарным значениям пункта назначения поездки.Это делается итеративно до тех пор, пока итоговые значения строк и столбцов не будут соответствовать удовлетворительному уровню точности.

# Setting origin (row) and destination (column) totals
TripOrig<-matrix(c(100,90,60,40,110))
TripDest<-matrix(c(60,110,50,100,80))

# Produce cost (impedance) matrix
CostMatrix<-matrix((rnorm(5*5,1,1.5))^2,5,5)

# Produce empty matrix for population with numbers of trips
TripInitial <- matrix(data=NA, nrow=5, ncol=5)

# Loop to populate initial trip matrix based on function

for(j in 1:5){
  for(i in 1:5){
    TripInitial[i,j] <-CostMatrix[i,j]^-0.9*TripOrig[i,1]*TripDest[j,1]
  }
}

# Set max iterations and the 'seed' trip matrix
MaxIter<-100
IterPrev<-TripInitial

# Bi-proportional fitting - iterate until convergence
for(i in 1:MaxIter) {

  EndIter<-i

  IterOrig<-matrix(data=NA, nrow=5, ncol=5)
  IterDest<-matrix(data=NA, nrow=5, ncol=5)

  for(j in 1:5){
    for(i in 1:5){
      IterOrig[i,j] = TripOrig[i,1]*IterPrev[i,j]/sum(IterPrev[i,])
    }
  }

  for(j in 1:5){
    for(i in 1:5){
      IterDest[i,j] = TripDest[j,1]*IterOrig[i,j]/sum(IterOrig[,j])
    }
  }

  RowTotals<-matrix(rowSums(IterDest))
  ColTotals<-matrix(colSums(IterDest))

  DiffOrig<-RowTotals-TripOrig
  DiffDest<-ColTotals-TripDest

  ConVal<-max(abs(DiffOrig),abs(DiffDest))

  IterPrev<-IterDest

  if (ConVal<.0001) {
    break
  }

}

Мой вопрос: могут ли эти необходимые матричные операции быть более эффективными при использовании семейства функций apply,как итерация по строкам и столбцам по отдельности не кажется вычислительно эффективной для больших матриц?