The expected value of the second pivot in gauss jordan eliminaton

Question

Say I have a matrix

x1 x2
x3 x4

With x1, x2, x3 and x4 randomly and uniformly drawn from the interval [0,1]

After I do gauss-jordan elimination, what is the expected value of the absolute value of the second pivot (the one on the x4 position)?

So far I seem to be getting a surprinsing answer: infinity

For the second pivot can be calculated as $x4-(x3/x1)*x2$. Therefore, the sought expected value is $E(x4-(x3/x1)*x2) = E(x4)-E(x2)*E(x3)*E(1/x1) = 0.5-0.5*0.5*E(1/x1)$.

But $E(1/x1)=\infty$ (see Expectation of 1/x, x uniform from 0 to 1).

Is there an error somewhere in my thinking? Where can I read a bit more about this?

Answer 1

Your computation appears to be correct to me. It's also supported by the following little bit of matlab code:

function gjtest(n)
mats = rand(4, n);

cells = zeros(n, 1);
for i = 1:n
    s = reshape(mats(:,i), 2, 2);
    s(2,:) = s(2,1)/s(1,1) * s(2,:);
    cells(i) = s(2,2);
end
mean(cells)
clf;
plot(cells); 
figure(gcf);

which produces, as mean values:

n       mean
1000     2.17...
10000    2.14...
100000   2.94...
1000000  3.100...
10000000 4.441...

which is not growing very fast, I admit, but certainly suggests that the mean is not $1$.

So you have a random variable with a very high mean, but the odd characteristic that its "population mean" tends to be a great deal smaller. I think this indicates lots of "skew" in the data (I apologize if I'm misusing statistics terms...it's been a long time). The plots from running those code fragments certainly suggest a very biased distribution. If you replace "plot" with "histogram" you can see what I mean.

In fact, let me add a little more. The pivot $X$ is (up to small perturbation), roughly $1/x_1$. How is this distributed? Well,

\begin{align} P(a-u < X < a+u) &= P(a-u < \frac{1}{x_1} < a + u)\\ &= P(\frac{1}{a+u} < x_1 < \frac{1}{a - u}\\ &= \frac{1}{a-u} - \frac{1}{a + u}\\ &= \frac{(a+u) - (a-u)}{(a-u)(a+u)}\\ &= \frac{2u}{a^2-u^2}\\ &\approx \frac{2u}{a^2}\\ \end{align} for small values of $u$. The pdf, gotten by taking a limit of difference quotients, therefore gives $$ p(x) = \frac{1}{a^2} $$ This explains why there are so few large values of $X$ -- the probability falls off quadratically with the target values. But it also explains why the expected value is infinite: if you look at $x p(x)$, whose integral is the mean of $X$, you find you're integrating $\frac{1}{x}$, which gives you a $\log$, and the integral diverges. Then again, it diverges about as slowly as something possibly can [I'm speaking informally here!], so it's no surprise that the values I was getting were in the $2$ to $5$ range. :)