Solve a matrix equation by Jacobi's method (Revised)

Question

Mathematical Background

Let A be an N by N matrix of real numbers, b a vector of N real numbers and x a vector N unknown real numbers. A matrix equation is Ax = b.

Jacobi's method is as follows: decompose A = D + R, where D is the matrix of diagonals, and R is the remaining entries.

if you make an initial guess solution x0, an improved solution is x1 = inverse(D) * (b - Rx) where all multiplications are matrix-vector multiplication and inverse(D) is the matrix inverse.

---

Problem Specification

• Input: Your complete program should accept as input the following data: the matrix A, the vector b, an initial guess x0, and an 'error' number e.
• Output: The program must output the minimum number of iterations such that the latest solution differs by the true solution, by at most e. This means each component of the vectors in absolute magnitude differ by at most e. You must use Jacobi's method for iterations.

How the data is inputted is your choice; it could be your own syntax on a command line, you could take input from a file, whatever you choose.

How the data is outputted is your choice; it can be written to a file, displayed in the command line, written as ASCII art, anything, as long as it is readable by a human.

Further Details

You are not given the true solution: how you calculate the true solution is up to you entirely. You may solve it by Cramer's rule for example, or computing an inverse directly. What matters is that you have a true solution to be able to compare to iterations.

Precision is an issue; some people's 'exact solutions' for comparison may differ. For the purposes of this code golf the exact solution must be true to 10 decimal places.

To be absolutely clear, if even one component of your present iteration solution exceeds its corresponding component in the true solution by e, then you need to keep iterating.

The upper limit to N varies based on what hardware you're using and how much time you're willing to spend running the program. For the purposes of this code golf, assume maximum N = 50.

Preconditions

When your program is called, you are free to assume that the following holds at all times:

• N > 1 and N < 51, i.e. you will never be given a scalar equation, always a matrix equation.
• All inputs are over the field of real numbers, and will never be complex.
• The matrix A is always such that the method converges to the true solution, such that you can always find a number of iterations to minimise the error (as defined above) below or equal to e.
• A is never the zero matrix or the identity matrix, i.e. there is one solution.

---

Test Cases

A = ((9, -2), (1, 3)), b = (3,4), x0 = (1,1), e = 0.04

The true solution is (0.586, 1.138). The first iteration gives x1 = (5/9, 1), differing by more than 0.04 from the true solution, by at least one component. Taking another iteration we find, x2 = (0.555, 1.148) which differs by less than 0.04 from (0.586, 1.138). Thus the output is

2

---

A = ((2, 3), (1, 4)), b = (2, -1), x0 = (2.7, -0.7), e = 1.0

In this case the true solution is (2.2, -0.8) and the initial guess x0 already has error less than e = 1.0, thus we output 0. That is, whenever you do not need to make an iteration, you simply output

0

---

Submission Assessment

This is code golf, with all standard loopholes hereby disallowed. The shortest correct complete program (or function), i.e. lowest number of bytes wins. It is discouraged to use things like Mathematica which wrap up a lot of the necessary steps into one function, but use any language you want.

Answer 1

APL (Dyalog), 78 68 65 49 bytes

Exactly the type of problem APL was created for.

-3 thanks to Erik the Outgolfer. -11 thanks to ngn.

Anonymous infix function. Takes A b e as left argument and x as right argument. Prints result to STDOUT as vertical unary using 1 as tally marks, followed by 0 as punctuation. This means that even a 0-result can be seen, being no 1s before the 0.

{⎕←∨/e<|⍵-b⌹⊃A b e←⍺:⍺∇D+.×b-⍵+.×⍨A-⌹D←⌹A×=/¨⍳⍴A}

Try it online!

Explanation in reading order

Notice how the code reads very similarly to the problem specification:

{…} on the given A, b, and e, and and the given x,
 ⎕← print
 ∨/ whether there is any truth in the statement that
 e< e is smaller than
 |⍵- the absolute value of x minus
 b⌹ b matrix-divided by
 ⊃A b e the first of A, b, and e (i.e. A)
 ←⍺ which are the left argument
 : and if so,
  ⍺∇ recurse on
  D+.× D matrix-times
  b- b minus
  ⍵+.×⍨ x, matrix multiplied by
  A- A minus
  ⌹D the inverse of D (the remaining entries)
  ← where D is
  A× A where
  =/¨ there are equal
  ⍳ coordinates for
  ⍴A the shape of A (i.e. the diagonal)

Step-by-step explanation

The actual order of execution right-to-left:

{…} anonymous function where ⍺ is A b e and ⍵ is x:
 A b c←⍺ split left argument into A, b, and e
 ⊃ pick the first (A)
 b⌹ matrix division with b (gives true value of x)
 ⍵- differences between current values of x and those
 | absolute values
 e< acceptable error less than those?
 ∨/ true for any? (lit. OR reduction)
 ⎕← print that Boolean to STDOUT
 : and if so:
  ⍴A shape of A
  ⍳ matrix of that shape where each cell is its own coordinates
  =/¨ for each cell, is the vertical and horizontal coordinates equal? (diagonal)
  A× multiply the cells of A with the that (extracts diagonal)
  ⌹ matrix inverse
  D← store in D (for Diagonal)
  ⌹ inverse (back to normal)
  A- subtract from A
  ⍵+.×⍨ matrix multiply (same thing as dot product, hence the .) that with x
  b- subtract that from b
  D+.× matrix product of D and that
  ⍺∇ apply this function with given A b e and that as new value of x

Answer 2

Python 3, 132 bytes

f=lambda A,b,x,e:e<l.norm(x-dot(l.inv(A),b))and 1+f(A,b,dot(l.inv(d(d(A))),b-dot(A-d(d(A)),x)),e)
from numpy import*
l=linalg
d=diag

Try it online!

Uses a recursive solution.

Answer 3

R, 138 bytes

function(A,x,b,e){R=A-(D=diag(diag(A)))
g=solve(A,b)
if(norm(t(x-g),"M")<e)T=0
while(norm((y=solve(D)%*%(b-R%*%x))-g,"M")>e){T=T+1;x=y}
T}

Try it online!

thanks to NikoNyrh for fixing a bug

It's also worth noting that there is an R package, Rlinsolve that contains a lsolve.jacobi function, returning a list with x (the solution) and iter (the number of iterations required), but I'm not sure that it does the correct computations.

Answer 4

Clojure, 212 198 196 bytes

#(let[E(range)I(iterate(fn[X](map(fn[r b d](/(- b(apply +(map * r X)))d))(map assoc % E(repeat 0))%2(map nth % E)))%3)](count(for[x I :while(not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))]x)))

Implemented without a matrix library, it iterates the process 1e9 times to get the correct answer. This wouldn't work on too ill-conditioned inputs but should work fine in practice.

Less golfed, I was quite happy with the expressions of R and D :) The first input % (A) has to be a vector, not a list so that assoc can be used.

(def f #(let[R(map assoc %(range)(repeat 0))
             D(map nth %(range))
             I(iterate(fn[X](map(fn[r x b d](/(- b(apply +(map * r x)))d))R(repeat X)%2 D))%3)]
          (->> I
               (take-while (fn[x](not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))))
               count)))