Rate Distorsion Theory with Time varying model

Question

I'm working with rate-distortion theory in a time-varying context, and I'm encountering a problem where the calculated rate goes below 0. I would like to get some insights into potential issues with my model and how to address them.

Here's a brief overview of the rate-distortion theory and the time-varying model I'm working with:

Rate-distortion theory is a fundamental framework in information theory that deals with the trade-off between the compression of a source and the resulting distortion after decompression. In a nutshell, the rate-distortion function quantifies the minimum rate (in bits per source symbol) needed to achieve a specific distortion level.

In a time-varying context, the underlying source probabilities and distortion levels may change over time. This adds complexity to the problem, as the rate-distortion function now depends on the evolving probabilities and distortion levels. The time-varying model I'm working with considers this dynamic aspect by incorporating time-varying source probabilities and distortion levels.

The main issue I'm encountering is that, under certain circumstances, the calculated rate goes below 0. This is an unexpected outcome, as the rate should be non-negative. I suspect there might be some issues or inconsistencies in my model or calculations that lead to this behavior.

Bellow is my python script which is somehow a representation of the Blahut-arimoto algorithm.

import numpy as np

def findOptimalChannel(s, time_periods=2):
    
    # Define source
    N = 5
    p_list = [
        [0.45, 0.4, 0.05, 0.05, 0.05],
        [0.35, 0.5, 0.05, 0.05, 0.05]
    ]

    # Initialize matrix
    A_list = []
    for t in range(time_periods):
        A = np.ones((N, N))
        for i in range(0, N):
            for j in range(i):
                A[i, j] = np.exp(s + t * 0.1)
                A[j, i] = np.exp(s + t * 0.1)
        A_list.append(A)

    # initial guess
    q = np.ones(N) / N

    Tu = 1
    Tl = 0
    max_iterations = 1000
    iterations = 0
    while (Tu - Tl > 0.001) and (iterations < max_iterations):
        iterations += 1
        for t in range(time_periods):
            A = A_list[t]
            p = p_list[t]
            b = A.dot(q)
            c = np.zeros(N)
            for k in range(N):
                for j in range(N):
                    c[k] += A[j, k] * p[j] / b[j]

            q = q * c

            Tu = -np.sum(q * np.log(c))
            Tl = -np.max(np.log(c))
    print('out')
    # Compute channel
    Q = np.zeros((N, N))
    for k in range(N):
        for j in range(N):
            Q[k, j] = A[j, k] * q[k] / b[j]

    # Compute distortion
    D = 0
    for j in range(N):
        dummy = 0
        for k in range(N):
            dummy += Q[k, j] * (1 - (j == k))
        D += p[j] * dummy

    # Compute rate
    R = s * D - np.sum(p * np.log(b)) - np.sum(q * np.log(c))

    return R, D, A_list, Q

raw = []
for k in range(len(s)):
    raw.append(findOptimalChannel(s[k]))