# Giant Babylonian Matrices

### Introduction

The Babylonians had a clever method for finding square roots. To find the square root of a number N, you start with a guess, then refine it by repeatedly evaluating a line of computer code (or papyrus scroll equivalent):

guess = (guess + N/guess) / 2


If the guess is smaller than sqrt(N), N/guess must be larger than sqrt(N), so averaging these two terms should get you closer to the true answer. It turns out this method converges fairly quickly to the correct answer:

N = 7
guess = 1
guess = 4.0
guess = 2.875
guess = 2.65489
guess = 2.64577
guess = 2.64575
guess = 2.64575


### Extension

Does an equivalent algorithm exist for matrices? Throwing all mathematical rigor aside, if we wish to find the 'square root' of a square matrix M, we may again start with a guess, then refine it by repeatedly evaluating a line of code:

guess = (guess + inverse(guess)*M) / 2


where inverse() is the matrix inverse function, and * is matrix multiplication. Instead of 1, our initial guess is an identity matrix of the appropriate size.

However, convergence is no longer guaranteed - the algorithm diverges for most (but not all) inputs. If M produces a convergent series, it is a Babylonian matrix. Furthermore, the converged value is actually the real square root of the matrix (i.e., guess*guess = M).

### Example

Here's an example using a 2x2 matrix, showing convergence of the series:

M = [[7,2],[3,4]]
guess = [[1.,      0.      ], [0.,       1.     ]]  <-- 2x2 identity matrix
guess = [[4.,      1.      ], [1.5,      2.5    ]]
guess = [[2.85294, 0.558824], [0.838235, 2.01471]]
guess = [[2.60335, 0.449328], [0.673992, 1.92936]]
guess = [[2.58956, 0.443035], [0.664553, 1.92501]]
guess = [[2.58952, 0.443016], [0.664523, 1.92499]]
guess = [[2.58952, 0.443016], [0.664523, 1.92499]]


After only a few iterations, the series converges to within 1e-9. After a while, it grows unstable as rounding errors add up, but the initial stability is demonstrated, which is sufficient.

Your job is to find a square matrix M, populated by the integers 1 through 9, which, when given an initial guess of an identity matrix of the proper size, converges successfully. Matrices should be submitted as comma-delimited text (or attached text files, if they are too large), with either brackets or newlines separating subsequent rows.

### Rules

Even though the matrix M and starting guess use integers (mainly for display convenience), all arithmetic should be conducted using double-precision floating point. Integer-specific operators (such as 3/2 evaluating to 1) are prohibited.

Divide-by-zero errors, singular matrices, and over/underflows cannot be used to coerce a sequence into 'converging'. The values must be proper floating-point numbers at all times to be considered to have converged.

The algorithm is considered to have converged if an iteration produces no change larger than 1e-9, in any individual cell.

You may use external mathematics libraries for matrix inversion (no need to reinvent the wheel).

### Win Condition

The winner is the largest submitted valid Babylonian matrix made of integers 1-9, after one week. You must also submit the code you used to find the matrix.

• So basically, get the code finished and fast, then run it for the rest of the week so that we can get the biggest matrix? This is not a good winning condition. Perhaps [fastest-code] would be appropriate for this (but then make sure you give a time limit AND something to make it fair across platforms with different power). May 23, 2014 at 7:01
• Good point. Perhaps this was not the most well-designed puzzle. I shall do better next time. May 23, 2014 at 8:24
• I see it differently as a math challenge disguised as a coding challenge, with the unstated goal of finding the rule that determines what makes the iteration converge or not, and using it to construct large Babylonian matrices.
– xnor
May 23, 2014 at 8:39
• @steveverrill, 0 isn't in 1-9. May 23, 2014 at 9:17

We produce such a matrix of any dimension n. It is simply a matrix of all 1's except 2's on the diagonal.

r = range(n)
M = [[1 + (i==j) for i in r] for j in r]


Example:

n=5
M=[[2 1 1 1 1]
[1 2 1 1 1]
[1 1 2 1 1]
[1 1 1 2 1]
[1 1 1 1 2]]

#Within 10 steps, the procedure converges to the matrix's square root
[[ 1.28989795  0.28989795  0.28989795  0.28989795  0.28989795]
[ 0.28989795  1.28989795  0.28989795  0.28989795  0.28989795]
[ 0.28989795  0.28989795  1.28989795  0.28989795  0.28989795]
[ 0.28989795  0.28989795  0.28989795  1.28989795  0.28989795]
[ 0.28989795  0.28989795  0.28989795  0.28989795  1.28989795]]


In general, the procedure converges for matrices with positive real eigenvalues. This is because if we work in the basis of eigenvectors of the original matrix, it and all its iterates are diagonal matrices, and we can view the procedure as applying the Babylonian method to each individual diagonal entry. The starting matrix is the still the identity in this basis.

Since the Babylonian method converges only if the original value was positive, the matrix version converges when every eigenvalue is real and positive. I'm not sure what happens for complex eigenvalues.

I constructed the given matrix to have positive eigenvalues as follows. First, start with the matrix of all ones, which is rank-1, and so has a single eigenvalue (of n) and all others 0. Then, add the identity, which shifts the eigenvalues up by 1 to get (n+1,1,1,...,1). Iterating makes the eigenvalues converge to (sqrt(n+1),1,1,...,1). Note that the 1's are unchanged by each iteration.

• In case the problem intended that all numbers 1-9 be used, an arbitrary matrix of 1's except diagonal >1 suffices. Sketch of proof: This matrix is symmetric and every principal minor has positive determinant (using the determinant lemma with update vectors of all ones), so by Sylvester's criterion, it is positive definite, which implies all positive real eigenvalues.
– xnor
May 23, 2014 at 7:57
• I see I was a bit hasty with this puzzle. Would it get more interesting if the matrix was required to be asymmetric, or to use only integers 2-9? May 23, 2014 at 8:26
• For assymetric, I suspect I can make the diagonal entries high and change a 1 to a 2. For integers 2-9, just double my original matrix. Honestly, I think any problem of this type is liable to be attackable by math rather than computer search, so any solution is likely to extend to arbitrarily large matrices.
– xnor
May 23, 2014 at 8:27
• I guess the puzzle is just broken. Maybe I should throw in the towel, eh? May 23, 2014 at 8:41
• I enjoyed it as a linear algebra problem, but it probably won't work as a coding challenge.
– xnor
May 23, 2014 at 8:43

Here's a 10x10 Babylonian matrix as a first target:

[[9, 8, 6, 5, 1, 1, 6, 7, 1, 6],
[8, 7, 4, 9, 9, 5, 4, 8, 5, 9],
[7, 1, 8, 8, 6, 8, 6, 6, 4, 1],
[7, 8, 9, 9, 5, 1, 1, 6, 9, 9],
[6, 9, 3, 4, 1, 9, 1, 7, 4, 9],
[9, 1, 1, 1, 5, 4, 5, 6, 7, 1],
[8, 6, 8, 2, 4, 4, 6, 8, 2, 2],
[5, 4, 9, 3, 3, 5, 5, 6, 4, 1],
[2, 4, 4, 8, 6, 8, 3, 3, 7, 6],
[3, 9, 9, 7, 7, 9, 6, 5, 1, 8]]


which converges to

[[3.20594, 1.60102, 0.237288, -0.0745693, -1.07811, -0.255869, -0.0182269, 1.70256, 0.355871, 1.11176],
[0.62321, 1.19922, -0.514309, 2.2604, 1.96759, 1.41485, -0.066537, 1.37868, -0.613093, 1.35704],
[0.919805, 0.453479, 2.54444, 0.717525, -0.00569383, 1.54527, 0.131113, 0.916813, 0.443482, -0.0251371],
[0.163383, -0.597711, 0.751538, 3.69952, 2.37728, -0.578767, 1.24375, -0.0922065, 0.749781, 0.623908],
[2.27446, 2.7975, 0.171454, -1.6263, -1.12639, 1.22686, -2.69745, 2.36239, 1.6901, 2.1349],
[0.669566, -1.57916, -0.722141, 1.35884, 2.43496, 2.25887, 2.01205, 0.241215, 0.342665, -0.714146],
[0.50365, 0.472729, 1.37547, 0.32966, 1.36842, 0.0550458, 2.90591, 0.520055, -0.257277, -0.347814],
[0.205106, 0.286162, 1.89821, 0.141372, 0.521793, 0.0638244, 1.65691, 1.04034, 0.81816, -0.30955],
[0.0679229, 1.01126, 0.553701, 0.627696, -0.0400446, 1.21608, -0.397756, 0.35419, 2.38639, 1.08573],
[-0.200249, 2.31035, 1.89905, -0.124084, 0.0398479, 0.454703, 1.08552, -0.152166, 0.408594, 2.31063]]


and was found with this Mathematica code:

n = 10;
M = RandomInteger[{1, 9}, {n, n}];
guess = IdentityMatrix[n] // N;
result = NestList[(# + Inverse[#].M)/2 &, guess, 32];
ListPlot[Transpose[Flatten /@ result], Joined -> True]
M // MatrixForm