Заполнение симметричной матрицы и поиск выражения максимально быстрым способом

Question

У меня есть n x n симметричная матрица G, чей элемент (i,j)th равен g(h(i),h(j)), с g(i,j) = g(j,i) и g(i,i) является постоянным для всех i. Здесь g принимает значения в реальной строке. В моем случае, скажем, g - это ядро ​​Гаусса.

Я попытался ввести матрицу следующим образом.

h = 0.01
operator = function(x,y){
  return(dnorm((x-y)/h))
}

a = 1; b = 0.5; n = 20000
x = seq(1,n,1)
vec = a*x + b*x^2 #example of h(x)
G = outer(vec, vec, FUN = operator)

Но здесь я вычисляю все записи матрицы, что не является необходимым. Достаточно только нижней треугольной матрицы и только одного элемента диагонали. Что я могу сделать, чтобы реализовать это? (Я думаю, что использование ifelse делает код медленным.)

Затем я хочу сделать следующее для некоторых двух n x 1 векторов a и b

a = rnorm(n); b = rcauchy(n)
s = rowSums(G)
sum(((a/s) * (G %*% (b/a)))^2)

Я знаю, что использование многоядерных процессоров с parallel делает код быстрее. Но я не знаю, как использовать это в моей обстановке. Любые предложения будут высоко оценены.

N.B. В разделе комментариев есть несколько предложений. Это, конечно, делает код быстрее, и поэтому я действительно ценю это. Но я ищу способ сделать этот кусок кода еще быстрее, если это возможно.

Answer 1

Извините, но ваши ожидания кажутся совершенно нереальными.Чтобы вычислить нижнюю часть матрицы 20k * 20k, нам нужно вычислить (n * nn) / 2 = 199990000 значений.Но,

system.time(dnorm((x-y)/h)) # 12.47 

Это самое простое вычисление в моей системе для случайных векторов такой длины занимает ~ 12 сек.Но нам все еще нужно создать большую симметричную матрицу.Поэтому я не понимаю, как это можно вычислить за ~ 0,005 с (используя R).

Из моего текущего тестирования ваш подход кажется самым быстрым способом сделать это в R.