Seeking help to implement a RAM matrix for structural equation models

[lf-araujo]

I have been trying to generate a minimal working optimization of the Reticular Action Model from McArdle 1984, with the general formula F (I-A)^{-1} S (I-A)^{-T} F^T. Where F is a filter matrix, S holds the covariances, and A the regression paths. The objective is to work on this until I get a general structural equation model optimizer in Nim as an exercise for myself, and maybe a future little package.

I figured I could start from the Rosenbrock example in the docs, and with the help from open ai, the farthest I got was:

import ggplotnim, numericalnim, benchy, std / [math,  sequtils]

# Define matrices F, A, and S
proc createMatrices(): (Tensor[float], Tensor[float], Tensor[float]) =
  let F = newTensor[float](3, 2) # Matrix F (3x2)
  F[0, 0] = 1.0
  F[0, 1] = 0.0
  F[1, 0] = 0.0
  F[1, 1] = 1.0
  F[2, 0] = 1.0
  F[2, 1] = 1.0

  let A = newTensor[float](2, 2) # Matrix A (2x2)
  A[0, 0] = 0.5
  A[0, 1] = 0.2
  A[1, 0] = 0.3
  A[1, 1] = 0.6

  let S = newTensor[float](2, 2) # Matrix S (2x2)
  S[0, 0] = 1.0
  S[0, 1] = 0.5
  S[1, 0] = 0.5
  S[1, 1] = 1.0

  return (F, A, S)

proc ramPath(theta: Tensor[float], F: Tensor[float], A: Tensor[float], S: Tensor[float]): float =
  let I = eye[float](2) # Identity matrix
  let invA = inv(I - A)  # Calculate (I - A)^{-1}
  result = F.dot(invA).dot(S).dot(invA.transpose()).dot(F.transpose())
  return result[0, 0]  # Return the scalar result

proc ramGradient(theta: Tensor[float], F: Tensor[float], A: Tensor[float], S: Tensor[float]): Tensor[float] =
  # Just an example of simple gradient calculations; actual gradient calculation can vary greatly
  let grad = newTensor[float](2)
  # Dummy gradient calculation, just using the parameter
  grad[0] = 2 * (ramPath(theta, F, A, S) - 1) * theta[0]
  grad[1] = 2 * (ramPath(theta, F, A, S) - 1) * theta[1]
  return grad

# Main process
proc main() =
  let (F, A, S) = createMatrices()
  let theta0 = [0.0, 0.0].toTensor()  # Initial guess
  let solutionLbfgs = lbfgs(ramPath, theta0, analyticGradient=ramGradient, extra=[F, A, S])
  echo "LBFGS Solution: ", solutionLbfgs

main()

I am still very far from a solution, right now I just needed help on identifying the correct procedure to generate the identity matrix correctly to substitute this dummy code: let I = eyefloat # Identity matrix.

Sorry if this is not the correct place to ask for this kind of help.

[HugoGranstrom]

Hi! :D I don't think this is a wrong place to ask, but maybe the #science channel on Discord/Matrix would be better as I'm not sure how many get notifications on these discussions. But I can answer here 😄

right now I just needed help on identifying the correct procedure to generate the identity matrix correctly to substitute this dummy code: let I = eyefloat # Identity matrix.

I don't really understand the question here. There is an eye function in arraymancer: https://mratsim.github.io/Arraymancer/special_matrices.html#eye%2Cvarargs%5Bint%5D
Is it the size of the matrix you are unsure about? Or what specifially is it that you need help with?