Seeking help to implement a RAM matrix for structural equation models #58

lf-araujo · 2024-12-15T16:19:07Z

lf-araujo
Dec 15, 2024

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 · 2024-12-15T19:48:08Z

HugoGranstrom
Dec 15, 2024
Maintainer

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?

2 replies

lf-araujo Dec 31, 2024
Author

Thank you!

I went back to it, fixed some errors in other parts of the code and changed optimization to steep. However, I am still running into non-conformable matrices in the formula calculation. I can't identify where the error is being produced, or I cant understand Nim's output in this case.

The code now looks like:

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

#     /  0.00         0.00    0.00    0.00    0.00  \
#     |  1.00         0.00    0.00    0.00    0.00  |
# A = |  \lam..2=0.33 0.00    0.00    0.00    0.00  |
#     |  \lam..3=0.33 0.00    0.00    0.00    0.00  |
#     \  0.00         0.00    0.00    0.00    0.00  /
# 
# 
#     /  sigma_s=1.00 0.00    0.00    0.00    0.00  \
#     |  0.00         e1=1.00 0.00    0.00    0.00  |
# S = |  0.00         0.00    e2=1.00 0.00    0.00  |
#     |  0.00         0.00    0.00    e3=1.00 0.00  |
#     \  0.00         0.00    0.00    0.00    0.00  /
# 
# 
#     /  0            1       0       0       0     \
# F = |  0            0       1       0       0     |
#     \  0            0       0       1       0     /

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

  # Matrix A (5x5)
  var A = newTensor[float](5, 5)
  A[0, 0] = 0.0
  A[0, 1] = 0.0
  A[0, 2] = 0.0
  A[0, 3] = 0.0
  A[0, 4] = 0.0
  A[1, 0] = 1.0
  A[1, 1] = 0.0
  A[1, 2] = 0.0
  A[1, 3] = 0.0
  A[1, 4] = 0.0
  A[2, 0] = 0.33
  A[2, 1] = 0.0
  A[2, 2] = 0.0
  A[2, 3] = 0.0
  A[2, 4] = 0.0
  A[3, 0] = 0.33
  A[3, 1] = 0.0
  A[3, 2] = 0.0
  A[3, 3] = 0.0
  A[3, 4] = 0.0
  A[4, 0] = 0.0
  A[4, 1] = 0.0
  A[4, 2] = 0.0
  A[4, 3] = 0.0
  A[4, 4] = 0.0

  # Matrix S (5x5)
  var S = newTensor[float](5, 5)
  S[0, 0] = 1.0
  S[0, 1] = 0.0
  S[0, 2] = 0.0
  S[0, 3] = 0.0
  S[0, 4] = 0.0
  S[1, 0] = 0.0
  S[1, 1] = 1.0
  S[1, 2] = 0.0
  S[1, 3] = 0.0
  S[1, 4] = 0.0
  S[2, 0] = 0.0
  S[2, 1] = 0.0
  S[2, 2] = 1.0
  S[2, 3] = 0.0
  S[2, 4] = 0.0
  S[3, 0] = 0.0
  S[3, 1] = 0.0
  S[3, 2] = 0.0
  S[3, 3] = 1.0
  S[3, 4] = 0.0
  S[4, 0] = 0.0
  S[4, 1] = 0.0
  S[4, 2] = 0.0
  S[4, 3] = 0.0
  S[4, 4] = 0.0

  return (F, A, S)

proc ramPath(theta: Tensor[float], F: Tensor[float], A: Tensor[float], S: Tensor[float]): float =
  let I = eye[float](5) # 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  # Return the scalar result

# Main process
proc main() =
  let (F, A, S) = createMatrices()
  # create a 5x5 theta 0 to.Tensor
  let theta0 = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1].toTensor() # Initial guess
  let solutionSteep = solutionSteepest(ramPath, theta0 )
  echo "Steep Solution: ", solutionSteep

main()

and the error:

/home/luis/.nimble/bin/nim c -r --verbosity\:0 --hint\[Processing\]\:off --excessiveStackTrace\:on /home/luis/Documents/OneDrive/Coding/Nim/Learning/MLestimator.nim
/home/luis/.nimble/pkgs2/ginger-0.6.1-01ae18b65875ed684d47d4b8f45be125d3621018/ginger.nim(3516, 24) Warning: Vega backend does not exist. [User]
/home/luis/.nimble/pkgs2/nimblas-0.3.1-e1ecdea4bb8176f12d66efd4aa0c7b3bea970027/nimblas/private/common.nim(52, 7) Hint: Using BLAS library matching pattern: lib(blas|cblas|openblas).so(|.3|.2|.1|.0) [User]
/home/luis/.nimble/pkgs2/nimlapack-0.3.1-fcb25795c6fb43f9251b7f34c3ad84c39e645afd/nimlapack.nim(20, 7) Hint: Using LAPACK library matching pattern: liblapack.so(|.3|.2|.1|.0) [User]
/home/luis/Documents/OneDrive/Coding/Nim/Learning/MLestimator.nim(101, 14) template/generic instantiation of `eye` from here
/home/luis/.nimble/pkgs2/arraymancer-0.7.33-c37fb5d98fce8661ade439d89d0a6a3c9d65c8fc/arraymancer/linear_algebra/special_matrices.nim(249, 9) template/generic instantiation of `set_diagonal` from here
/home/luis/.nimble/pkgs2/arraymancer-0.7.33-c37fb5d98fce8661ade439d89d0a6a3c9d65c8fc/arraymancer/linear_algebra/special_matrices.nim(172, 26) Error: type mismatch
Expression: &"Diagonal index ({k=}) exceeds input matrix height ({a.shape[0]})"
  [1] "Diagonal index ({k=}) exceeds input matrix height ({a.shape[0]})": string

Expected one of (first mismatch at [position]):
[1] func `&`(a, b: DynamicStackArray): DynamicStackArray
[1] func `&`[T](a: DynamicStackArray[T]; value: T): DynamicStackArray[T]
[1] proc `&`(x, y: char): string
[1] proc `&`(x: char; y: string): string
[1] proc `&`[T](x, y: sink seq[T]): seq[T]
[1] proc `&`[T](x: sink seq[T]; y: sink T): seq[T]
[2] proc `&`(x, y: string): string
[2] proc `&`(x: string; y: char): string
[2] proc `&`[T](x: sink T; y: sink seq[T]): seq[T]

I can correct the matrix sizes, but I don't think I am passing incorrect ones. What am I missing?

HugoGranstrom Dec 31, 2024
Maintainer

Looking at the eye call it seems like it expects two inputs. Not sure ehy you don't get an error for that though :/ It also seems like there is an identity function for generating the identity matrix: https://github.com/mratsim/arraymancer/blob/master/src%2Farraymancer%2Flinear_algebra%2Fspecial_matrices.nim#L228
Try using that instead and see if it works.

lf-araujo · 2024-12-31T18:41:25Z

lf-araujo
Dec 31, 2024
Author

Thanks. It really helped. It seems I made some progress, however the optimization algo procedure (steepestDescend()) seems to require a float and my formula returns a Tensor. I will look into how to optimize based on the covariance matrix (result). Also, this only works with openblas in my machine nim c -d:blas=openblas -r MLestimator.nim

import ggplotnim, numericalnim, arraymancer, benchy

#     /  0.00         0.00    0.00    0.00    0.00  \
#     |  1.00         0.00    0.00    0.00    0.00  |
# A = |  \lam..2=0.33 0.00    0.00    0.00    0.00  |
#     |  \lam..3=0.33 0.00    0.00    0.00    0.00  |
#     \  0.00         0.00    0.00    0.00    0.00  /
# 
# 
#     /  sigma_s=1.00 0.00    0.00    0.00    0.00  \
#     |  0.00         e1=1.00 0.00    0.00    0.00  |
# S = |  0.00         0.00    e2=1.00 0.00    0.00  |
#     |  0.00         0.00    0.00    e3=1.00 0.00  |
#     \  0.00         0.00    0.00    0.00    0.00  /
# 
# 
#     /  0            1       0       0       0     \
# F = |  0            0       1       0       0     |
#     \  0            0       0       1       0     /
# Model Covariance matrix = F (I-A)^{-1} S (I-A)^{-T} F^T


proc ramPath(F: Tensor[float], A: Tensor[float], S: Tensor[float]): Tensor[float] =
  var F = [[0.0, 1.0, 0.0, 0.0, 0.0],
           [0.0, 0.0, 1.0, 0.0, 0.0],
           [0.0, 0.0, 0.0, 1.0, 0.0]].toTensor().reshape(3, 5) # Reshape to 3x5

  # Matrix A (5x5)
  var A = [[0.0, 0.0, 0.0, 0.0, 0.0],
           [1.0, 0.0, 0.0, 0.0, 0.0],
           [0.33, 0.0, 0.0, 0.0, 0.0],
           [0.33, 0.0, 0.0, 0.0, 0.0],
           [0.0, 0.0, 0.0, 0.0, 0.0]].toTensor()
  # Matrix S (5x5)
  var S = [[1.0, 0.0, 0.0, 0.0, 0.0],
           [0.0, 1.0, 0.0, 0.0, 0.0],
           [0.0, 0.0, 1.0, 0.0, 0.0],
           [0.0, 0.0, 0.0, 1.0, 0.0],
           [0.0, 0.0, 0.0, 0.0, 0.0]].toTensor() 
  let I = [[1.0, 0.0, 0.0, 0.0, 0.0],
           [0.0, 1.0, 0.0, 0.0, 0.0],
           [0.0, 0.0, 1.0, 0.0, 0.0],
           [0.0, 0.0, 0.0, 1.0, 0.0],
           [0.0, 0.0, 0.0, 0.0, 1.0]].toTensor()
  let invA = pinv(I - A)  # Calculate (I - A)^{-1}
  var result = F * invA * S * invA.transpose() * F.transpose()
  return result

let theta0 = [
   [1,1,1,1,1],
        [1,1,1,1,1],
        [1,1,1,1,1],
        [1,1,1,1,1],
        [1,1,1,1,1]].toTensor()

let solutionSteepest = steepestDescent(ramPath, theta0)
echo "Steepest: ", solutionSteepest

/home/luis/Documents/OneDrive/Coding/Nim/Learning/MLestimator.nim(56, 39) Error: type mismatch
Expression: steepestDescent(ramPath, theta0)
  [1] ramPath: proc (F: Tensor[system.float], A: Tensor[system.float], S: Tensor[system.float]): Tensor[system.float]{.gcsafe.}
  [2] theta0: Tensor[system.int]

Expected one of (first mismatch at [position]):
[1] proc steepestDescent[U; T: not Tensor](f: proc (x: Tensor[U]): T; x0: Tensor[U];
    options: OptimOptions[U, StandardOptions] = steepestDescentOptions[U]();
    analyticGradient: proc (x: Tensor[U]): Tensor[T] = nil): Tensor[U]
[1] proc steepest_descent(deriv: proc (x: float64): float64; start: float64;
                      gamma: float64 = 0.01; precision: float64 = 0.00001;
                      max_iters: Natural = 1000): float64

5 replies

HugoGranstrom Dec 31, 2024
Maintainer

numericalnim expects a function of the form f(Tensor[float]): Tensor[float]. Where the Tensors are 1-dimensional Tensors. So you have to rewrite the function to not take in F, A & S but it should take just one input.
I'll look up if we have a good documentation on this and link it

HugoGranstrom Dec 31, 2024
Maintainer

https://scinim.github.io/getting-started/numerical_methods/optimization.html

See this tutorial and read the linked Further reading documentation of numericalnim

lf-araujo Dec 31, 2024
Author

Thank you. This is it for an initial attempt. The code below seems to be optimizing the covariance matrix, and basically returning approx what I have given as starting values, which is what I expected. My next steps will be to provide observed data and optimize from the expected covariance matrix to the observed covariance matrix. Thanks for the explanations.

import ggplotnim, numericalnim, arraymancer, benchy

proc ramPath(theta: Tensor[float]): float =
  # Reshape theta to separate F, A, and S matrices
  var F = theta[0..14].reshape(3, 5)
  var A = theta[15..39].reshape(5, 5)
  var S = theta[40..64].reshape(5, 5)

  let I = [[1.0, 0.0, 0.0, 0.0, 0.0],
                         [0.0, 1.0, 0.0, 0.0, 0.0],
                         [0.0, 0.0, 1.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0, 1.0, 0.0],
                         [0.0, 0.0, 0.0, 0.0, 1.0]].toTensor()

  let invA = pinv(I - A)  # Calculate (I - A)^{-1}
  var result = F * invA * S * invA.transpose() * F.transpose()
  #echo "Result: ", result
  return result[0,0]


let theta0 = [
  # Initial values for F
  0.0, 1.0, 0.0, 0.0, 0.0,
  0.0, 0.0, 1.0, 0.0, 0.0,
  0.0, 0.0, 0.0, 1.0, 0.0,

  # Initial values for A
  0.0, 0.0, 0.0, 0.0, 0.0,
  1.0, 0.0, 0.0, 0.0, 0.0,
  0.33, 0.0, 0.0, 0.0, 0.0,
  0.33, 0.0, 0.0, 0.0, 0.0,
  0.0, 0.0, 0.0, 0.0, 0.0,

  # Initial values for S
  1.0, 0.0, 0.0, 0.0, 0.0,
  0.0, 1.0, 0.0, 0.0, 0.0,
  0.0, 0.0, 1.0, 0.0, 0.0,
  0.0, 0.0, 0.0, 1.0, 0.0,
  0.0, 0.0, 0.0, 0.0, 0.0].toTensor()

let solutionSteepest = steepestDescent(ramPath, theta0)
echo "Steepest: ", solutionSteepest

Hint: /home/luis/Documents/OneDrive/Coding/Nim/Learning/MLestimator [Exec]
Steepest: Tensor[system.float] of shape "[65]" on backend "Cpu"
    -0.00162307     0.00186136    0.000185275    0.000185275              0              0              0              1              0              0              0              0              0              1              0     -0.0845598      -0.226571     -0.0173847     -0.0173847              0       0.903068      -0.396114     -0.0139398     -0.0139398              0       0.336354      0.0232686    0.000237471    0.000237471              0       0.336354      0.0232686    0.000237471    0.000237471              0              0              0              0              0              0       0.936159     -0.0982497     0.00652027     0.00652027              0     -0.0982497       0.779634      0.0129567      0.0129567              0     0.00652027      0.0129567       0.998819    -0.00118125              0     0.00652027      0.0129567    -0.00118125       0.998819              0              0              0              0              0              0

HugoGranstrom Jan 1, 2025
Maintainer

Glad to have helped, hope it works as expected :D. Now that you have got it working using steepestdescent it should work with the other optimization algorithms as well 😃 lbfgs is the best one in most cases

HugoGranstrom Jan 1, 2025
Maintainer

Just ping me if you get into more problem ^^

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

SciNim

Seeking help to implement a RAM matrix for structural equation models #58

{{title}}

Replies: 2 comments 7 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

SciNim

Seeking help to implement a RAM matrix for structural equation models #58

lf-araujo Dec 15, 2024

Replies: 2 comments · 7 replies

HugoGranstrom Dec 15, 2024 Maintainer

lf-araujo Dec 31, 2024 Author

HugoGranstrom Dec 31, 2024 Maintainer

lf-araujo Dec 31, 2024 Author

HugoGranstrom Dec 31, 2024 Maintainer

HugoGranstrom Dec 31, 2024 Maintainer

lf-araujo Dec 31, 2024 Author

HugoGranstrom Jan 1, 2025 Maintainer

HugoGranstrom Jan 1, 2025 Maintainer

lf-araujo
Dec 15, 2024

Replies: 2 comments 7 replies

HugoGranstrom
Dec 15, 2024
Maintainer

lf-araujo Dec 31, 2024
Author

HugoGranstrom Dec 31, 2024
Maintainer

lf-araujo
Dec 31, 2024
Author

HugoGranstrom Dec 31, 2024
Maintainer

HugoGranstrom Dec 31, 2024
Maintainer

lf-araujo Dec 31, 2024
Author

HugoGranstrom Jan 1, 2025
Maintainer

HugoGranstrom Jan 1, 2025
Maintainer