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())