Return to Answer

% Helper functions
    % Binary entropy
fn_h = @(p) -p.*log2(p) - (1-p).*log2(1-p); 
    % MI across mtx_bc when v_distn is input
fn_I = @(mtx_bc,v_distn) fn_h(v_distn(1)) + fn_h(v_distn*mtx_bc(:,1)) ...
    - sum(sum(-log2(diag(v_distn)*mtx_bc).*(diag(v_distn)*mtx_bc)));
    % Channel matrix when P(out=0|in=0)=pa, P(out=0|in=1)=pb
fn_mtxBC = @(pa,pb) [pa, 1-pa; pb, 1-pb];

% Set params
d_r_L = 0.4; 
d_r_R = 0.4;
d_xp = 0.146102; % solution to 1-H(p) = 0.4
mtxBSC = fn_mtxBC(d_xp, 1-d_xp);

% Search 
while true
    mtxL = fn_mtxBC(rand, rand);
    mtxR = fn_mtxBC(rand, rand);
    v_d = (mtxL'\[0.5, 0.5]')';
    if (abs(sum(v_d)-1) > 0.001 || ...
        min(v_d) < 0)
        continue
    end
    if(fn_I(mtxL, v_d)      > 0.4 || ...
       fn_I(mtxR, v_d*mtxL) > 0.4)
        continue;
    end
    fprintf('+\n');
    if fn_I(mtxL*mtxR, v_d) > fn_I(mtxBSC*mtxBSC, [0.5, 0.5])
        break;
    end
end

% Helper functions
    % Binary entropy
fn_h = @(p) -p.*log2(p) - (1-p).*log2(1-p); 
    % MI across mtx_bc when v_distn is input
fn_I = @(mtx_bc,v_distn) fn_h(v_distn(1)) + fn_h(v_distn*mtx_bc(:,1)) ...
    - nansum(nansum(-log2(diag(v_distn)*mtx_bc).*(diag(v_distn)*mtx_bc)));
    % Channel matrix when P(out=0|in=0)=pa, P(out=0|in=1)=pb
fn_mtxBC = @(pa,pb) [pa, 1-pa; pb, 1-pb];

% Set params
d_r_L = 0.4; 
d_r_R = 0.4;
d_xp = 0.146102; % solution to 1-H(p) = 0.4
mtxBSC = fn_mtxBC(d_xp, 1-d_xp);

% Search 
while true
    mtxL = fn_mtxBC(rand, rand);
    mtxR = fn_mtxBC(rand, rand);
    v_d = (mtxL'\[0.5, 0.5]')';
    if (abs(sum(v_d)-1) > 0.001 || ...
        min(v_d) < 0)
        continue
    end
    if(fn_I(mtxL, v_d)      > 0.4 || ...
       fn_I(mtxR, v_d*mtxL) > 0.4)
        continue;
    end
    fprintf('+\n');
    if fn_I(mtxL*mtxR, v_d) > fn_I(mtxBSC*mtxBSC, [0.5, 0.5])
        break;
    end
end

An easier-to-investigate (and arguably more interesting) problem is one identical to yours but omitting the fourth constraint, that $X\sim B(1/2)$. But I could not find a satisfying answer for this one either.

An easier-to-investigate (and arguably more interesting) problem is one identical to yours that omits the fourth constraint that $X\sim B(1/2)$. But I could not find an easy path towards an answer for this either.

No. You can get a higher $I(U;V)$ using asymmetric channels. Below I construct a counterexample, but first a more succinct restatement of the question.

% Helper functions
    % Binary entropy
fn_h = @(p) -p.*log2(p) - (1-p).*log2(1-p); 
    % MI across mtx_bc when v_distn is input
fn_I = @(mtx_bc,v_distn) fn_h(v_distn(1)) + fn_h(v_distn*mtx_bc(:,1)) ...
    - sum(sum(-log2(diag(v_distn)*mtx_bc).*(diag(v_distn)*mtx_bc)));
    % Channel matrix when P(out=0|in=0)=pa, P(out=0|in=1)=pb
fn_mtxBC = @(pa,pb) (pa>=pb)*[pa, 1-pa; pb, 1-pb]; fn_p = @(mtx) mtx(:,1);

% Set params
d_r_L = 0.4; 
d_r_R = 0.4;
d_xp = 0.146102; % solution to 1-H(p) = 0.4
mtxBSC = fn_mtxBC(d_xp, 1-d_xp);

% Search 
while true
    mtxL = fn_mtxBC(rand, rand);
    mtxR = fn_mtxBC(rand, rand);
    v_d = (mtxL'\[0.5, 0.5]')';
    if (abs(sum(v_d)-1) > 0.001 || ...
        min(v_d) < 0)
        continue
    end
    if(fn_I(mtxL, v_d)      > 0.4 || ...
       fn_I(mtxR, v_d*mtxL) > 0.4)
        continue;
    end
    fprintf('+\n');
    if fn_I(mtxL*mtxR, v_d) > fn_I(mtxBSC*mtxBSC, [0.5, 0.5])
        break;
    end
end

No. You can get a higher $I(U;V)$ using asymmetric channels. Below I construct a counterexample, but first a more succinct restatement of the question.

% Helper functions
    % Binary entropy
fn_h = @(p) -p.*log2(p) - (1-p).*log2(1-p); 
    % MI across mtx_bc when v_distn is input
fn_I = @(mtx_bc,v_distn) fn_h(v_distn(1)) + fn_h(v_distn*mtx_bc(:,1)) ...
    - sum(sum(-log2(diag(v_distn)*mtx_bc).*(diag(v_distn)*mtx_bc)));
    % Channel matrix when P(out=0|in=0)=pa, P(out=0|in=1)=pb
fn_mtxBC = @(pa,pb) [pa, 1-pa; pb, 1-pb];

% Set params
d_r_L = 0.4; 
d_r_R = 0.4;
d_xp = 0.146102; % solution to 1-H(p) = 0.4
mtxBSC = fn_mtxBC(d_xp, 1-d_xp);

% Search 
while true
    mtxL = fn_mtxBC(rand, rand);
    mtxR = fn_mtxBC(rand, rand);
    v_d = (mtxL'\[0.5, 0.5]')';
    if (abs(sum(v_d)-1) > 0.001 || ...
        min(v_d) < 0)
        continue
    end
    if(fn_I(mtxL, v_d)      > 0.4 || ...
       fn_I(mtxR, v_d*mtxL) > 0.4)
        continue;
    end
    fprintf('+\n');
    if fn_I(mtxL*mtxR, v_d) > fn_I(mtxBSC*mtxBSC, [0.5, 0.5])
        break;
    end
end