Inspired by this.
If we start with a matrix [[1, 1], [1, 1]] and we perform iterative proportional fitting with real-valued non-zero growth factors [a, b], I've found that the equilibrium matrix has the form:
[2a^2 / (a + b), 2ab / (a + b)]
[2ab / (a + b), 2b^2 / (a + b)]
By growth factors, I mean the desired row and column sums are given by the original ones times those factors.
What if we have a different initial matrix?
If your initial matrix is all 5s instead of all 1s, then all the entries in the equilibrium matrix are 5 times larger.
If you scale the initial matrix columns by positive integers (m, n), then in the equilibrium matrix, the columns are scaled by (m, n).
If you scale the initial rows by constants, you don't get anything like that though. Which is weird, because the iterative-proportional-fitting procedure is pretty symmetric between rows and columns. I guess in my code I take turns scaling rows (first) and columns (second) over and over, so that's one source of asymmetry.
Anyway, let's say we have an initial matrix [[m, m], [n, n]]. What's the equilibrium matrix for different growth factors (a, b)? Weirdly it only converges when the growth factors are equal, {a = b}. Otherwise, for example, the top left matrix entry alternates between two small rational values. When the growth factors are equal, the equilibrium matrix is the initial matrix scaled by the growth factor, i.e. [[am, am], [an, an]].
Okay, we've seen fully-constant initial matrices, and matrices with constant columns, and matrices with constant rows. Now for something spicier. Let's investigate a matrix that has all 1s except for the top left entry. Let's have a small change of notation: [a, b] no longer refer to growth factors. Instead our initial matrix is [[a, b], [c, d]], so that we can say {b=c=d=1}.
I've figured out symbolic expressions for the equilibrium matrix when we have growth factors [1/p, 1] for {p} a positive integer. I don't know of any reason why these formulas wouldn't work with {p} drawn from a more general classes of numbers, but I figured out the formulas using {p} over positive integers, and that's the only place where I've checked the formulas, and that's the only place where I'll assert their correctness now.
All the entries can be written as solutions of quadratic equations with leading terms {px^2} and rational coefficients. For given values of {a} and {p}, all the quadratic equations associated with the four entries of the equilibrium matrix have the same discriminant.
The top left entry of the equilibrium matrix is
((p + a^2 + a/2 - 1/2) - sqrt(D)) / (p(a - 1))
for discriminant
D = p^2 + (2a^2 + a - 1)p + ((a+1)/2)^2
The anti-diagonal of the equilibrium matrix is constant. Both values are given by
(-(p + (a+1)/2) + sqrt(D)) / (p(a -1))
The bottom right entry of the equilibrium matrix is
((2ap + a/2 - p + 1/2) - sqrt(D)) / (p(a - 1))
My method for figuring out these formulas mostly amounted to finding symbolic forms for entries of equilibrium matrices with defined {a} and {p}, and gradually noticing linear and quadratic patterns in constants, generalizing repeatedly from special cases. All pattern recognition, no proofs. But the formulas are correct. And if you can explain why, that'd be cool to see.
Next I'd like to try finding symbolic expressions for the equilibrium matrix terms when we have independent growth rates, or perhaps when we have non-unitary values for both {a} and {d} in the initial matrix, but those sound pretty hard, and I'm just going to savor the accomplishment of finding the above formulas for a while.
Oh! I can figure out growth rates [1, 1/p] instead of [1/p, 1]. That sounds like a manageable challenge. Please hold!
...
For growth rates [1, 1/p] and initial matrix [[a, 1], [1, 1]], IPF gives us an equilibrium matrix with top left entry
((a^2 + a/2 - 1/2)p + 1 - sqrt(D)) / ((a - 1)p)
for discriminant
D = ((a + 1)/2)^2 p^2 + (2a^2 + a - 1) p + 1
The antidiagonal is constant, and both values are given by
(-1 * (((a + 1)/2)p + 1) + sqrt(D)) / ((a - 1)p)
and the bottom right entry of the equilibrium matrix is
(p(a + 1)/2 + (2a - 1) - sqrt(D)) / ((a - 1)p)
Nice.
I just checked, and these formulas still work if your second growth rate isn't the reciprocal of an integer. You can use a growth rates of [1, 6/1] or [1, 5/8] or whatever. Just set {p} equal to the reciprocal of your second growth rate and they'll still predict correctly. And the case is likewise the set of formulas before these, for growth rates [1/p, 1]. You can use positive reals in the first slot, just set {p} to the reciprocal of the first growth rate.
...
I want to take small steps to generalize from here. I'm thinking that I could either look at initial matrices of the form [[a, 1], [1, a]] or I could look at growth rates {c * [1/p, 1]} and {c * [1, 1/p]}. And then eventually I'll get to independent diagonal values and independent growth rates. Or maybe that's already equivalent to independent growth rates? Because two real numbers will be related by a multiplicative constant {c}, namely their ratio. And {p} will be the inverse of one of them. I think that might be everything. Welp.
Let's start with an initial matrix of the form [[a, 1], [1, a]] and growth factors [1/p, 1], for {p} a positive integer. I fully expect that the formulas will continue to generalize to {1/p} being any positive real number.
For a=3, the discriminant term of entries in the equilibrium matrices has the form
D = p^2 + 34p + 1
The top left entry is given by
(p + 17 - sqrt(D)) / 4p
The anti-diagonal is constant with both values given by
(-p - 1 + sqrt(D)) / 4p
And the bottom right entry is given by
(17p + 1 - sqrt(D)) / 4p
.
For a = 5, my first guess was that the discriminant was
D = 1/4 p^2 + 49/2 p + 1/4
but that looks a little funny, and you can move a constant between the discriminant and the rest of the expression. So let's also consider four times the previous expression:
D = p^2 + 98p + 1
If we use that 4x discriminant, then the upper left entry of the equilibrium matrix has the form
(p + 49 - sqrt(D)) / 8p
which isn't so bad. Let's stick with this discriminant. The anti-diagonal is constant, with values of the form
(-p + -1 + sqrt(D)) / 8p
and the bottom right entry looks like
(49p + 1 - sqrt(D)) / 8p
Smashing. Let's do {a = 7} and then we'll guess and check the full form for general {a} and then move on to something more challenging.
For {a = 7}, the discriminant term in the equilibrium matrix entries is
D = p^2 + 194p + 1
The top left entry is
(p + 97 - sqrt(D)) / 12p
the anti-diagonal is
(-p - 1 + sqrt(D)) / 12p
and the bottom right entry is
(97p + 1 - sqrt(D)) / 12p
.
So, generalizing over a=(3, 5, 7), it looks like the general discriminant term of the equilibrium matrix entries generated by IPF from a starting matrix [[a, 1], [1, a]] and growth factors [1/p, 1] is:
p^2 + (4a^2 - 2)p + 1
We'll check that of course, but I've heard of crazier things.
The top left entry generalized to depend on {a} looks like
(p + (2a^2 - 1) - sqrt(D)) / (2(a - 1)p)
The anti-diagonal is
(-p - 1 + sqrt(D)) / (2(a - 1)p)
and the bottom right entry is
((2a^2 - 1)p + 1 - sqrt(D)) / (2(a - 1)p)
...
I had the thought that if I do analyze initial matrix
[[a, 1], [1, a]]
and the
[[a, 1], [1, 2a]]
and then
[[a, 1], [1, 3a]]
and maybe a few more, then I could try to generalize over those formulas and hopefully get a formula that's free in the coefficient in front of {a} on the bottom left, and then if that generalizes past the positive integers, then I've basically got independent upper-left and bottom-right entries. And then hopefully the simple patterns continue where you can multiply columns or rows by constants. And then we've got the whole damn thing solved. For one growth factor at least.
I'm starting to think that this problem is too boring to solve.
...
Maybe I can power through it. Just do 5 minutes a day.
Let's work on the case of initial matrix [[a, 1], [1, 2a]], growth rates [1/p, 1].
For a = 3, the discriminant term is
D = ?
The top left entry is
(p + ? - sqrt(D)) / (?p)
The antidiagonal is constant with value
(-(p + 1) + sqrt(D)) / (?p)
The bottom right entry is
(?p + ? - sqrt(D)) / (?p)
...
No comments:
Post a Comment