Search code examples
rsum

How to write the denominator summation?


I first sample a random matrix and its eigenvalues and eigenvectors.

n <- 2000

#Sample GOE random matrix
A <- matrix(rnorm(n*n, mean=0, sd=1), n, n) 
G <- (A + t(A))/sqrt(2*n)
ev <- eigen(G)
l <- ev$values
v <- ev$vectors

and I also generate random vectors uniformly distributed on the sphere

#size of multivariate distribution
mean <- rep(0, n) 
var <- diag(n)

#make this example reproducible
set.seed(101)

#simulate bivariate normal distribution
initial <- mvrnorm(n=10, mu=mean, Sigma=var)

#normalized the first possible initial value, the initial data uniformly distributed on the sphere
x_0 <- initial[1, ]/norm(initial[1, ], type="2") 

So x_0 is a random vector uniformly distributed on the sphere and each column of v are eigenvector of the matrix.

Now I want to define the following function

![enter image description here

My code is as follows. But my question is how to write the denominator.

h10 <- x_0 %*% v[, n]

sum((x_0 %*% v[, i])^2) #How to fix this one?

h1 <- function(t) {
  abs(h10)*exp(-2*l[n]*t)/(sqrt())
}

Solution

  • Here is a way. Break the computations into numer and denom and use R's vectorized ability to compute the matrix products, power and square root.
    The function is scalar on t. So I also vectorize (unrelated to the above) the scalar function on "t".

    #n <- 2000
    n <- 20
    
    #Sample GOE random matrix
    A <- matrix(rnorm(n*n, mean=0, sd=1), n, n) 
    G <- (A + t(A))/sqrt(2*n)
    ev <- eigen(G)
    l <- ev$values
    v <- ev$vectors
    
    #size of multivariate distribution
    mean <- rep(0, n) 
    var <- diag(n)
    
    #make this example reproducible
    set.seed(101)
    
    #simulate bivariate normal distribution
    initial <- MASS::mvrnorm(n=10, mu=mean, Sigma=var)
    
    #normalized the first possible initial value, the initial data uniformly distributed on the sphere
    x_0 <- initial[1, ]/norm(initial[1, ], type="2")
    
    h1_scalar <- function(t, x0, v, lambda) {
      h10 <- x0 %*% v[, n]
      numer <- abs(h10) * exp(-2*lambda[n] * t)
      h <- sum((x0 %*% v)^2)
      denom <- h * exp(-4*lambda * t)
      c(numer)/sqrt(denom)
    }
    h1 <- Vectorize(h1_scalar, "t")
    
    h1(1, x_0, v, l)
    #>              [,1]
    #>  [1,] 13.62478196
    #>  [2,] 12.53497558
    #>  [3,]  7.50725898
    #>  [4,]  5.66588613
    #>  [5,]  4.22446682
    #>  [6,]  2.71595400
    #>  [7,]  1.97947406
    #>  [8,]  1.13463313
    #>  [9,]  0.91133245
    #> [10,]  0.75639112
    #> [11,]  0.57187518
    #> [12,]  0.36569767
    #> [13,]  0.29825876
    #> [14,]  0.20393953
    #> [15,]  0.18553709
    #> [16,]  0.12584079
    #> [17,]  0.08098017
    #> [18,]  0.04032765
    #> [19,]  0.03189976
    #> [20,]  0.01920207
    
    h1(1:10, x_0, v, l)
    #>              [,1]         [,2]         [,3]         [,4]         [,5]
    #>  [1,] 13.62478196 9.667431e+03 6.859502e+06 4.867143e+09 3.453470e+12
    #>  [2,] 12.53497558 8.182744e+03 5.341637e+06 3.486983e+09 2.276278e+12
    #>  [3,]  7.50725898 2.935045e+03 1.147488e+06 4.486229e+08 1.753940e+11
    #>  [4,]  5.66588613 1.671813e+03 4.932958e+05 1.455550e+08 4.294841e+10
    #>  [5,]  4.22446682 9.293852e+02 2.044653e+05 4.498249e+07 9.896175e+09
    #>  [6,]  2.71595400 3.841464e+02 5.433393e+04 7.685029e+06 1.086976e+09
    #>  [7,]  1.97947406 2.040570e+02 2.103553e+04 2.168479e+06 2.235409e+08
    #>  [8,]  1.13463313 6.704446e+01 3.961597e+03 2.340872e+05 1.383201e+07
    #>  [9,]  0.91133245 4.325194e+01 2.052742e+03 9.742338e+04 4.623725e+06
    #> [10,]  0.75639112 2.979510e+01 1.173662e+03 4.623188e+04 1.821126e+06
    #> [11,]  0.57187518 1.703156e+01 5.072332e+02 1.510640e+04 4.498980e+05
    #> [12,]  0.36569767 6.964603e+00 1.326388e+02 2.526066e+03 4.810817e+04
    #> [13,]  0.29825876 4.632745e+00 7.195874e+01 1.117709e+03 1.736097e+04
    #> [14,]  0.20393953 2.165982e+00 2.300425e+01 2.443214e+02 2.594866e+03
    #> [15,]  0.18553709 1.792724e+00 1.732192e+01 1.673705e+02 1.617192e+03
    #> [16,]  0.12584079 8.246978e-01 5.404658e+00 3.541943e+01 2.321213e+02
    #> [17,]  0.08098017 3.415146e-01 1.440257e+00 6.073940e+00 2.561540e+01
    #> [18,]  0.04032765 8.469502e-02 1.778741e-01 3.735663e-01 7.845536e-01
    #> [19,]  0.03189976 5.299400e-02 8.803717e-02 1.462532e-01 2.429656e-01
    #> [20,]  0.01920207 1.920207e-02 1.920207e-02 1.920207e-02 1.920207e-02
    #>               [,6]         [,7]         [,8]         [,9]        [,10]
    #>  [1,] 2.450401e+15 1.738676e+18 1.233673e+21 8.753500e+23 6.211024e+26
    #>  [2,] 1.485938e+15 9.700100e+17 6.332157e+20 4.133587e+23 2.698377e+26
    #>  [3,] 6.857222e+13 2.680906e+16 1.048129e+19 4.097776e+21 1.602070e+24
    #>  [4,] 1.267263e+13 3.739268e+15 1.103333e+18 3.255564e+20 9.606077e+22
    #>  [5,] 2.177164e+12 4.789775e+14 1.053753e+17 2.318264e+19 5.100194e+21
    #>  [6,] 1.537426e+11 2.174546e+13 3.075693e+15 4.350281e+17 6.153068e+19
    #>  [7,] 2.304404e+10 2.375530e+12 2.448850e+14 2.524434e+16 2.602350e+18
    #>  [8,] 8.173208e+08 4.829476e+10 2.853694e+12 1.686222e+14 9.963737e+15
    #>  [9,] 2.194425e+08 1.041477e+10 4.942861e+11 2.345887e+13 1.113361e+15
    #> [10,] 7.173619e+07 2.825769e+09 1.113102e+11 4.384635e+12 1.727157e+14
    #> [11,] 1.339884e+07 3.990437e+08 1.188430e+10 3.539378e+11 1.054096e+13
    #> [12,] 9.162056e+05 1.744886e+07 3.323084e+08 6.328713e+09 1.205284e+11
    #> [13,] 2.696616e+05 4.188555e+06 6.505930e+07 1.010542e+09 1.569639e+10
    #> [14,] 2.755930e+04 2.926992e+05 3.108672e+06 3.301629e+07 3.506563e+08
    #> [15,] 1.562587e+04 1.509826e+05 1.458846e+06 1.409588e+07 1.361993e+08
    #> [16,] 1.521207e+03 9.969234e+03 6.533339e+04 4.281625e+05 2.805963e+06
    #> [17,] 1.080269e+02 4.555776e+02 1.921290e+03 8.102585e+03 3.417073e+04
    #> [18,] 1.647698e+00 3.460449e+00 7.267539e+00 1.526308e+01 3.205510e+01
    #> [19,] 4.036306e-01 6.705381e-01 1.113943e+00 1.850556e+00 3.074266e+00
    #> [20,] 1.920207e-02 1.920207e-02 1.920207e-02 1.920207e-02 1.920207e-02
    

    Created on 2022-11-28 with reprex v2.0.2