We will try to adjust the adaptation of the noise-based with a smoother choice of Vnoise. For the correlation, we know that:

Vtask : 1 / sqrt(1 / Vnoise^2 * (1 / rmach - 1) - 1);
Vmach : 1 / sqrt(1 / Vnoise^2 * (1 / rtask - 1) - 1);

We would like the coefficient of variation of all the costs in the matrix to be Vmax:

V1 : Vtask;
V2 : Vmach;
V3 : Vnoise;
Vtotal : sqrt(V1^2*V2^2*V3^2 + V1^2*V2^2 + V1^2*V3^2 + V2^2*V3^2 + V1^2 + V2^2 + V3^2);

Let’s solve this with maxima:

assume(rtask >= 0, rmach >= 0, rtask <= 1, rmach <= 1);
assume(Vnoise > 0, Vmax > 0);
declare(rtask, constant);
declare(rmach, constant);
declare(Vmax, constant);
r : solve([Vtotal^2 - Vmax^2], [Vnoise]);
factor(r[4]);

Tests show that this gives positive values for the CV:

Vmax : 1;
rtask : 0.5;
rmach : 0.1;
r : solve([Vtotal^2 - Vmax^2], [Vnoise]);
Vnoise : r[4];
Vtask : 1 / sqrt(1 / Vnoise^2 * (1 / rmach - 1) - 1);
Vmach : 1 / sqrt(1 / Vnoise^2 * (1 / rtask - 1) - 1);
float([Vtask, Vmach, Vnoise]);

We found that the following choice of Vnoise was good:

Vmax <- 1
rtask <- 0.5
rmach <- 0.1
N1 <- -sqrt((rtask-rmach)^2*Vmax^4+
              2*(rtask*(rmach-1)^2+rmach*(rtask-1)^2)*Vmax^2+(rmach*rtask-1)^2)
N2 <- -(2*rtask*rmach-rtask-rmach)*Vmax^2-rmach*rtask+1
D <- 2*rtask*rmach*(Vmax^2+1)
(Vnoise <- sqrt((N1+N2)/D))
## [1] 0.563283
(Vtask <- 1 / sqrt(1 / Vnoise^2 * (1 / rmach - 1) - 1))
## [1] 0.1911608
(Vmach <- 1 / sqrt(1 / Vnoise^2 * (1 / rtask - 1) - 1))
## [1] 0.6817227

Let’s how Vtask and Vmach differ to the previous way:

k <- Vmax / sqrt(Vmax^2 + 1)
(Vnoise <- min(Vmax, sqrt(1 / max(rtask, rmach) - 1) * k))
## [1] 0.7071068
(Vtask <- Vnoise / sqrt(1 / rmach - Vnoise^2 - 1))
## [1] 0.2425356
(Vmach <- Vnoise / sqrt(1 / rtask - Vnoise^2 - 1))
## [1] 1

This is close. By symmetry, we now focus on Vtask (and also Vnoise):

Vmax <- 1
variations <- NULL
for (rtask in seq(1e-3, 1-1e-3, 0.03)) {
  for (rmach in seq(1e-3, 1-1e-3, 0.03)) {
    N1 <- -sqrt((rtask-rmach)^2*Vmax^4+
                  2*(rtask*(rmach-1)^2+rmach*(rtask-1)^2)*Vmax^2+(rmach*rtask-1)^2)
    N2 <- -(2*rtask*rmach-rtask-rmach)*Vmax^2-rmach*rtask+1
    D <- 2*rtask*rmach*(Vmax^2+1)
    Vnoise <- sqrt((N1+N2)/D)
    Vtask <- 1 / sqrt(1 / Vnoise^2 * (1 / rmach - 1) - 1)
    variations <- rbind(variations, data.frame(rtask=rtask, rmach=rmach,
                                               Vnoise=Vnoise, Vtask=Vtask))
  }
}

For Vnoise:

library(ggplot2)
ggplot(variations, aes(x=rtask, y=rmach, fill=Vnoise, z=Vnoise)) +
  geom_tile() +
  geom_contour() +
  scale_fill_continuous(limits=c(0, 1))

For Vtask:

ggplot(variations, aes(x=rtask, y=rmach, fill=Vtask, z=Vtask)) +
  geom_tile() +
  geom_contour() +
  scale_fill_continuous(limits=c(0, 1))

Now with the previous method:

Vmax <- 1
variations <- NULL
for (rtask in seq(1e-3, 1-1e-3, 0.03)) {
  for (rmach in seq(1e-3, 1-1e-3, 0.03)) {
    k <- Vmax / sqrt(Vmax^2 + 1)
    Vnoise <- min(Vmax, sqrt(1 / max(rtask, rmach) - 1) * k)
    Vtask <- Vnoise / sqrt(1 / rmach - Vnoise^2 - 1)
    variations <- rbind(variations, data.frame(rtask=rtask, rmach=rmach,
                                               Vnoise=Vnoise, Vtask=Vtask))
  }
}

For Vnoise:

library(ggplot2)
ggplot(variations, aes(x=rtask, y=rmach, fill=Vnoise, z=Vnoise)) +
  geom_tile() +
  geom_contour() +
  scale_fill_continuous(limits=c(0, 1))

For Vtask:

ggplot(variations, aes(x=rtask, y=rmach, fill=Vtask, z=Vtask)) +
  geom_tile() +
  geom_contour() +
  scale_fill_continuous(limits=c(0, 1))

The result is smoother with an equivalent objective: limit the overall coefficient of variation.