Thank you for the data-driven perspective. I hadn’t thought of that approach.
I modified your snippet to show the mean absolute value of the gradient too.
import torch
import torch.nn.functional as F
import sys

a = torch.randn(1000,1000, requires_grad=True)
b = a
print (f"in: {a.std().item():.4f}")
for i in range(100):
    l = torch.nn.Linear(1000,1000, bias=False)
    torch.nn.init.xavier_normal_(l.weight, float(sys.argv[2]))
    b = getattr(F, sys.argv[1])(l(b))
    if i % 10 == 0:
        print (f"out: {b.std().item():.4f}", end=" ")
        a.grad = None
        b.sum().backward(retain_graph=True)
        print (f"grad: {a.grad.abs().mean().item():.4f}")

A few results


sigmoid seems stable with any gain > 0, but the gradients vanish pretty fast.
The more layers you have the higher the gain you will need.


tanh seems stable with pretty much any gain > 1
With gain 5/3 the output stabilises at ~.65, but the gradients start to explode after around 10 layers
Gain 1.1 works much better, giving output std stable around 0.30 and grads that are much more stable though they do grow slowly


softsign with gain 1 has slowly vanishing output and gradients
Gain > 1 reduces the vanishing, but higher values eventually cause the gradients to explode
The higher the gains the faster the gradients explode as you add layers


relu seems to be inherently less stable than the others, but it works OK with gain ~= sqrt(2)


selu only works with gain 1 and gives output std ~= 1.00, but the grad slowly explodes after 10-20 layers