Problems in calculating the log absolute Jacobian determinant of the CorrCholeskyFactor transformation

[Junyi-Que]

The way TransformVariables.jl constructs the Cholesky factor of a correlation seems similar to what's given here https://mc-stan.org/docs/reference-manual/correlation-matrix-transform.html, with the only difference being $z = tanh(y/2)$ instead of $z = tanh(y)$. The log absolute Jacobian determinant should be $$\sum{i=1}^{K-1}\sum{j=i+1}^{K}\frac{K-i-1}{2}\ln[1-(z{i,j})^2]+\sum{i=1}^{K-1}\sum{j=i+1}^{K}\ln\frac{\partial z{i,j}}{\partial y{i,j}}$$ If $K=3$, it's $$\frac{1}{2}(\ln[1-(z{1,2})^2]+\ln[1-(z{1,3})^2])+\sum{i=1}^{2}\sum{j=i+1}^{3}\ln\frac{\partial z{i,j}}{\partial y_{i,j}}$$

But the calculation in the library misses $\frac{1}{2}\ln[1-(z_{1,2})^2]$ when K=3, and this problem generalizes to larger Ks. The code below illustrates this problem for $K=3$.

t = as(Array,CorrCholeskyFactor(3))
x =  [-0.12833806457027003, 0.24721188811516634, -0.24578748441593382]
logjac_pack = TransformVariables.transform_and_logjac(t,x)[2]

#= 
Calculate the change-of-variable adjustment for LKJ transformation
=#
function den_adj(V::Vector{T}) where T
    N = length(V);
    K_X = Int(sqrt(2*N+1/4)+1/2);
    J_X = 0
    count = 1
    @inbounds for j in 2:K_X
        for i in 1:j-1
            yij =V[count]
            J_X += (K_X-i-1)/2*(log(4)+yij-2*log(exp(yij)+1))+log(2)+yij-2*log(1+exp(yij))
            count += 1
        end
    end
    return J_X
end

logjac_man = den_adj(x)

isequal(logjac_pack,den_adj(x)-1/2* (log(4)+x[1]-2*log(exp(x[1])+1)))  # this returns true

[tpapp]

Thanks for looking into this. I thought the code accounted for the difference, and we also test for this. Currently the tests use ForwardDiff, but we should probably move to FD. See eg

using TransformVariables, FiniteDifferences, LinearAlgebra
t = CorrCholeskyFactor(3)
x =  [-0.12833806457027003, 0.24721188811516634, -0.24578748441593382]
_, lj = TransformVariables.transform_and_logjac(t, x)
above_diag(M::AbstractMatrix) = [M[i, j] for (i, j) in map(Tuple, eachindex(M)) if i < j]
j_fd = jacobian(central_fdm(5, 1), x -> above_diag(transform(t, x)), x)[1]
isapprox(logabsdet(j_fd)[1], lj)

tests true.

[Junyi-Que]

Dear Tamas,

Sorry for the false alarm. I think what I computed was the log absolute Jacobian determinant of transforming to the correlation matrix instead of the Cholesky factor. Thank you for your patience.

Best, Junyi