Negative z_vec in Recursive Logit Model

[JunzeYang]

Hi Matt,
I'm currently using the recursive logit model to calculate the probabilities of transitioning from upstream to downstream links given a destination. However, I encountered an issue while solving for $z$ in the equation $z=Mz+b$ using _compute_exp_value_function. Specifically, I found that $z$ can sometimes take on negative values, which is unexpected since $z$ represents the exponential of the value function and should theoretically be greater than zero.

For example, in the Sioux Falls network, if I set all $v(a\mid k)$ to $-0.5$ and choose Node $3$ as the destination, there will be negative values in z_vec.

Do you have any suggestions on how to address this issue?

Thank you!

[m-richards]

Hi @JunzeYang, do you have a code example you could share? It's been quite some time since I thought deeply about this problem / last looked at the code and that might help me look at this a bit more quickly. It does not surprise me though that cases such as this are possible, I ran into situations where $(I-M)z=b$ was non invertible. I seem to recall running into non real solutions as well. My memory is pretty hazy at this point, I think in Tien mai's code there may have been a numerical capping of negative values, but I never found a mathematical justification to do so.

[JunzeYang]

Thank you very much for your prompt reply. Below is the code I use to calculate $z$:

def cal_exp_value_func(M:sparse.csr_matrix) -> numpy.ndarray:
    nrows = M.shape[0]
    rhs = sparse.lil_matrix((nrows, 1))
    rhs[-1, 0] = 1
    a_mat = sparse.identity(nrows) - M
    z = splinalg.spsolve(a_mat, rhs)
    realmin = np.finfo(np.float64).tiny
    z[z < realmin] = realmin
    z_abs = np.abs(z)
    z_abs = z_abs.reshape(nrows, 1)
    gap = (rhs - a_mat.dot(z_abs)).sum()
    print(f"gap of z: {gap}\n")
    return z_abs

I suspect this issue might be related to the built-in matrix inversion algorithm in Scipy. When I set large values to $v(a\mid k)$, the negative values in $z$ approach zero. However, when $v(a\mid k)$ is relatively small, the minimum negative values in $z$ may be far from zero, such as -6.

I have reviewed Mai's code, and it seems that they arbitrarily set a lower bound and change negative elements in $z$ to this bound without detailed discussion on this issue. Perhaps one feasible approach is to ensure $v(a\mid k)$ is large enough?

Below is Mai's code on computing $z$. Thank you again.

function [expV, boolV0] = getExpV(M)
    boolV0 = 1;
    [n m]= size(M);
    b = zeros(n,1);
    b(n)=1;
    b_sp = sparse(b);    
    I = speye(size(M));  
    A = I - M; 
    Z = A\b_sp;
    
    % test "values smaller than 0"
    minele = min(Z);
    if minele == 0 || minele < OptimizeConstant.NUM_ERROR
        boolV0 = 0;
    end  
    
    % set values smaller than real min to realmin
    Z(Z<realmin) = realmin;
    
    % set negative values (if any) to positive ones
    Zabs = abs(Z); 
    
    %test "still solution of system after making neg elements positive"
    resNorm = norm(A * Zabs - b_sp);
    if resNorm > OptimizeConstant.RESIDUAL
        boolV0 = 0;
    end

    expV = full(Zabs);
end

[m-richards]

Perhaps one feasible approach is to ensure $v(a\mid k)$ is large enough?

Given that utility is a relative concept and only differences in utility matter, if this is enough such that the solution $z$ is non-negative then that seems fine - to me it certainly more theoretically sound than the matlab snippet if it works. Reading the snippet (thank's for pulling this out, would have taken me quite a while to find it again), it seems like they intially tried arbitrarily taking the positive value and then superceded this with the realmin condition (if the realmin condition is applied, then Z should be strictly positive) as a workaround.