title: Filling up a symmetric matrix and finding an expression in a fastest possible way
body: <p>I have a <code>n x n</code> symmetric matrix <code>G</code>, whose <code>(i,j)th</code> element is <code>g(h(i),h(j))</code>, with <code>g(i,j) = g(j,i)</code> and <code>g(i,i)</code> is constant for all <code>i</code>. Here <code>g</code> takes values in real line. In my case, say, <code>g</code> is the Gaussian kernel. </p>
  
  <p>I have tried to input the matrix as follows.</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>But, here I am calculating all the entries of the matrix, which is not necessary. Only lower triangular matrix and only one element of the diagonal should suffice. What can I do to implement that? (I think using <code>ifelse</code> makes the code slow.)</p>
  
  <p>Then I want to do the following for some two <code>n x 1</code> vectors <code>a</code> and <code>b</code> </p>
  
  <pre><code>a = rnorm(n); b = rcauchy(n)
  s = rowSums(G)
  sum(((a/s) * (G %*% (b/a)))^2)
  
  </code></pre>
  
  <p>I know that using multi-cores with <code>parallel</code> makes code faster. But I do not know how to use that in my setting. Any suggestions will be much appreciated. </p>
  
  <p>N.B. There are some suggestions in the comment section. Those, of course, makes the code faster and thus I really appreciate those. But I am looking for a way which can make this whole piece of code even faster, if possible. </p>