Generating orthogonal latent variables

Don’t quite understand the context of your problem but do you want something like this (in this case they are orthonormal but could just make them orthogonal too)?

import pymc as pm 
import numpy as np

with pm.Model():
  
  vec = pm.MvNormal("mv", 0, [[1,0],[0,1]], size=2)
  
  X = vec[0,:]
  Y = vec[1,:]
  normX = pm.math.sqrt((X**2).sum())
  
  Z = X/normX
  W = Y-(Y*Z).sum()*Z
  normW = pm.math.sqrt((W**2).sum())
  W = W/normW

  print((W*W).sum().eval())
  print((Z*Z).sum().eval())
  print((W*Z).sum().eval())