Eliminating rounding errors from matrix

Question

My problem is that I have a Matrix where the sum of all rows, and the sum of all columns is zero. All numbers are rounded to x decimals.

Then I multiply the entire matrix with a number between 0 and 1 (eg. 1/6) and round all the numbers to x decimals. Now I cannot be sure that the sum of the rows and columns will be zero. I want the sums to be zero again with the smallest possible adjustment (or at least very small adjustment)

Does there exist a algoritm that can fix such a problem?

Example (very simple): matrix:

    200  -200  0

    400  400  -800

   -600 -200  800

round2( (1/6)*matrix)

33.33  -33.33  0   

66.67  66.67   -133.33

-100   -33.33  133.33

Answer 1

This is not a solution; just a more mathematical description of what you are trying to achieve (without judging if this is the right thing to do):

Since you are rounding all the numbers to x decimals, we can treat these numbers as integers (just multiply them by 10^x).

Now, you are trying to solve the following problem:

Given the matrix

A11+Adj11   A12+Adj12   ...   A1n+Adj1n
A21+Adj21   A22+Adj22   ...   A2n+Adj2n
A31+Adj31   A32+Adj32   ...   A3n+Adj3n
...         ...         ...   ...
Am1+Adjm1   Am2+Adjm2   ...   Amn+Adjmn

Where A11..Amn are constant integers,

Find integers Adj11...Adjmn

Minimizing sum(abs(Adjxy))

(or maybe you prefer: Minimizing sum((Adjxy)^2)

Subject to:

- for each row m: Adjm1+Adjm2+...+Adjmn = - (Am1+Am2+...+Amn)
- for each col n: Adj1n+Adj2n+...+Adjmn = - (A1n+A2n+...+Amn)

This is an integer programming problem, with m*n variables and m+n constrains. The function that you are trying to minimize is not linear.

I'm afraid that this problem is far from trivial. I believe that you should better post it on https://math.stackexchange.com/

Answer 2

What you're experiencing here is essentially a precision error. There's nothing you can do unless you don't round at all. This is similar to saving a photo as a 256 color image. You're losing information (essentially precision; due to discretization) and there's nothing you can do. For pictures, there are algorithms to make images appear smoother/closer to the original (e.g. dithering), but you don't have such things for single value numbers.

Possible solutions (actually just one with two different ways to visualize):

• Only round for display. It should be possible for the user to interpret that numbers are truncated/rounded. In your example it should be obvious that 6.67 would actually be 6.66666....

• Don't round at all and just truncate the numbers after a fixed number of decimals (add ... if needed; this is actually similar to the other solution).

In general, if you'd like to solve linear equations (or math in general), always use the available (and reasonable; performance wise) datatype with the biggest precision available, usually being single or double precision values. Otherwise you're introducing error margins getting worse and worse the more you calculate with them.

Answer 3

You can't eliminate rounding when you're working with floating point numbers. Your best solution might be to stick with integers in the matrix, then apply the final 1/6 to the result.

Answer 4

A common way to ensure that small rounding errors will not lead to a big error on the sum is to check that the error does not become too large at each partial sum.

With a one dimension vector [a[1], a[2], ..., a[n]], you can compute the partial sums [a[1], a[1]+a[2], ..., a[1]+a[2]+...+a[n]], multiply it and then restore the good vector by substracting the previous cell to the current one: [a[1]*b, (a[1]+a[2])*b-a[1]*b, ..., (a[1]+a[2]+...+a[n])*b-(a[1]+a[2]+...+a[n-1])*b]. By using this trick, the error on any partial sum is no more than 10^(-x).

You can adapt this method for a 2 dimensions matrix with the 3 following procedures:

partial_sum(M) =
  for i = 0 to n-1 do
    for j = 1 to m-1 do
      M[i][j] += M[i][j-1]
    done
  done
  for i = 0 to n-1 do
    for j = 1 to m-1 do
      M[j][i] += M[j-1][i]
    done
  done

multiply(M, a) =
  for i = 0 to n-1 do
    for j = 0 to m-1 do
      M[i][j] *= a
    done
  done

restore(M) =
  for i = 0 to n-1 do
    for j = 1 to m-1 do
      M[i][j] -= M[i][j-1]
    done
  done
  for i = 0 to n-1 do
    for j = 1 to m-1 do
      M[j][i] -= M[j-1][i]
    done
  done