Can We Solve Linear Algebra Problems at Extreme Scale and Low Precisions? | Hacker News Reader

Can We Solve Linear Algebra Problems at Extreme Scale and Low Precisions?

nhigham.com

95 pointsstncls25 comments

The article covers ways to obtain smaller error bounds, mainly for testing algorithm stability in practice.

The methods proposed are very varied in nature: blocked algorithms, x87 80-bit arithmetic, AVX2/ARMv8-A fused multiply-add (FMA) operations, and probabilistic bounds.

While the article is otherwise very good, I am unsure about one of those in 2021: 80-bit arithmetic. The corresponding x87 instructions are deprecated by Intel, and their performance is disastrous. For example, the complex algorithms implemented in software by glibc to compute transcendental functions (sin, cos, exp, ...) are faster and more accurate than their single instruction x87 FPU counterparts (FSIN, FCOS, ...)!

I would be very curious to see, in a few algorithms, a comparison of x87 arithmetic with double-doubles and even low-precision MPFR. MPFR, by the way, may have its place here if we're talking about sacrificing some performance for higher precision.

The trascendental instructions are slow and have low accuracy, but the "normal" x87 fp ops are still fast (except for the lack of vectorization).

You're right, I didn't know that. In a quick "long double" vs "double" microbenchmark, I see a negligible penalty for additions and comparisons, and 37% penalty when multiplications are involved.

Yeah. The problem with them is that using double-double based methods gives you 104 bits of precision that works with vectorization (and therefore is faster than x87).

I assume double-double would be the emulated __float128 type in GCC (or hand written equivalent)? Can GCC vectorize it or it needs to be done by hand?

__float128 is emulated quad-precision, I think.

Not exactly, the mantissa is increased, but the exponent is the same as standard double precision, so it is not an IEEE quad.

Agreed, I've seen better results from algorithmic approaches.

For example, Takahashi for bounded-error on extremely large dot product. Depending on use case even doing a pseudo-Takahashi where you only merge the error back once every N iterations can unblock instruction-level parallelism letting the performance hit be only around 2x relative to naive dot products, but with still effectively the same guarantees.

I think they were only recommending x87 for a final summation

Higham's writing is very readable and his book on numerical linear algebra is approachable and has information that's hard to find elsewhere.

His point about probabilistic error bounds for roundoff errors is excellent! The usual bounds are worst-case, and assume all roundoff errors complement each other. That doesn't happen in practice: the errors average out -- it's readily observable.

As Higham points out, this means intuitively that error magnitudes go from O(n) to O(sqrt(n)) which is huge for large n. I wasn't aware of his very recent joint work with Theo Mary that proves this under some conditions.

I've heard the phenomenon from your second paragraph called "the miracle of mitigating errors" in the context of old-school (pre-GPS) transoceanic navigation.

One interesting question is how to implement linear algebra for the case of binary variables x_i in {0,1} or x_i in {-1,0,+1} as opposed to having continuous variables. Of course, it will be somewhat different theory but the general properties should be analogous to those of standard linear algebra. For example, dot product should express the notion of orthogonality, there should be well defined distance etc. I am not sure that such a theory exists at all but if yes then I am curious if there exist implementations.

You might find what you're looking for if you google "linear systems over finite fields" or "integer linear systems".

Mixed integer programming?

Here is a browser version https://observablehq.com/@tomlarkworthy/mip

It exists. Linear algebra is defined over finite fields. In the case of {0, 1} that is Z_2 and the case of { -1, 0, + 1 } that is Z_3.

There are algorithms which fail over GF(2). Gram-schmidt, for instance.

Discrete mathematics theory is possibly more difficult than the set of complex numbers. Case in point: linear optimization across the reals/complex field is solved with the simplex algorithm.

Boolean optimization in contrast, is NP-complete since its equivalent to the knapsack problem (0 for "don't pack the object" and 1 for "do pack the object". Each object has a space it takes up + a value. Optimizing on value is NP-complete). Literally more difficult from a computational perspective than its Real/Complex analogue, and __provably__ so.

I'm enjoying my personal studies into Binary Decision Diagrams: the data-structures that are used in Verilog / VHDL synthesizers, libraries, CPU-design, verification (etc. etc.) to encode boolean truth tables and yes... tackle NP-complete and even #P-complete problems (#P complete is "determine the number of solutions" to an NP complete problem. Knights Tours aren't NP complete but... see the paper for an example of #P-like counting problems: "The Number of Knight's Tours Equals 33,439,123,484,294": https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.45...).

The boolean truth table for a 64x64-bit multiplication is 128-bits of input and 128-bits of output. That's 2^253 bytes of storage under a naive scheme, which is clearly unworkable. Binary Decision Diagrams make it possible to store, compute, and "think" about these functions at a high level. (Ex: counting solutions, finding the 1st solution, combining truth tables together, etc. etc.)

Knuth's TAOCP Volume 4A has been a good introduction to the theory, but Knuth's writing style is so difficult sometimes. I'm probably going to buy a supplemental textbook on the subject...

> Discrete mathematics theory is possibly more difficult than the set of complex numbers.

Probably it is so because in continuous case there is the luxury of having enough points between any other points. In discrete case you are much more constrained, for example, you cannot choose a point between 0 and 1, and what is the length of the diagonal of a binary cube or angle between its two hyper-planes? Sometimes these notions can be naturally defined but in other cases the formal theory is not so natural.

By the way, an interesting question is how complex boolean numbers could be defined (naturally).

> By the way, an interesting question is how complex boolean numbers could be defined (naturally).

Ever work with GF(2^x) extension fields?

For a good practical application of "binary complex numbers" (aka: GF(2^x) extension fields), see Reed Solomon error correction codes. EDIT: NASA's Reed Solomon code tutorial is an excellent, "casual", introduction to the subject. (Very little math theory involved: mostly sticks to the "final math" so to speak, focusing on the Reed Solomon error correction aspects). Once you understand how Reed Solomon codes are used, you can then more easily go into the theory of how it relates to GF extension fields, and how those GF-extension fields are similar to complex numbers made up of 0s and 1s.

Very interesting thank you. Is there an area of study that considers bit strings as one dimensional stochastic processes?

I'd have to imagine that's "just" a 2-valued discrete Fourier transform / discrete cosine transform.

>By the way, an interesting question is how complex boolean numbers could be defined (naturally).

Booleans act as a weird sort of field where there's no addition¹ and both ∧ and ∨ act as multiplication operators with F≡0 for ∧ and T≡0 for ∨. From this, you can build up the 3×3 tables for each relation to derive that if we add a third value to the boolean arithmetic, say, i, that

T ∨ i = T

F ∨ i = i

i ∨ i = F

and

F ∧ i = F

T ∧ i = i

i ∧ i = T

No other values are possible without destroying the group operations on ∧ and ∨ since we can't have a multiplication operation between two non-zero values give zero and within the non-zero members of the group, we cannot have repeated elements in any row or column.

———

1. This isn't strictly true since we can actually define two possible addition operators, each associated with ∧ and ∨. The associated addition for ∧ would be exclusive or (A + B = T iff A ≠ B) and the associated addition for ∨ would be if and only if (A + B = T iff A = B). The extension of these groups to include i is left as an exercise for the reader, but yes, there is only one possible solution.

You're trying to build up an extension field but you're just not quite getting there.

I'll try to describe the GF(2^2) field, aka, the GF(4) field, consisting of two binary numbers. In GF-world, we use the variable "x" instead of "i" to represent the "complex-like" arithmetic.

This means that every "complex boolean" (aka: GF(2^2) number) is a two-bit number. Lets name the bits a and b. That would be "number = a + bx", where "x" is the polynomial variable (aka: think of it "like i" from the complex world).

As such, we are actually dealing with a mapping of 2-bit numbers into 2-bit numbers. So our options have grossly increased over what you think.

------------

Because GF(4) is so tiny, I'll just give you a GF(4) algebra over these numbers and skip over all the theory. A GF(2^n) system can be constructed using a prime polynomial: all "numbers" are to be taken modulo the prime polynomial. For GF(4), one such field is constructed with the prime-number x^2 + x + 1 (aka: 111).

Yes, we have 2-bit numbers, but our modulus exists in the 3-bit space. Don't worry about it. Its needed for the math to work.

Addition is fortunately a simple XOR. So 00 + 11 == 11 (where + is the simple XOR operator).

Multiplication is... really complicated. We start with the "carryless multiplication" to be consistent with our XOR-based definition of add/subtract. Its basically described as "multiply, but throw away the carries". So 11 * 11 == 101 (the 2nd bit's carry is thrown away). But that's not enough: as per the theory of fields, we only reach a field if we... modulus it by a prime number (which is 111).

Its easier if I just gave you the table.

00 * anything == 0

01 * X == X

Okay, now we get complicated. All multiplications here on out "overflow", and need to be reduced modulo 111. This works because "addition" is defined as XOR, so we have the so called "Carryless multiply" operation as our primitive, to remain consistent with the addition operator. (aka: throw away carries, don't propagate them)

10 * 00 and 10 * 01 have been defined as earlier.

10 * 10 == 100 % 111 == 11

10 * 11 == 110 % 111 == 01

11 * 11 == 101 % 111 == 10

Given this multiplication table, we can see some interesting properties.

01 * 10 == 10

10 * 10 == 11

11 * 10 == 01

So 10^0 == "1", 10^1 == 10, 10^2 == 11, 10^3 == "1".

Which gives us the similar "looping" behavior of "i" in complex-number world. It also means that logarithms, square roots, cube-roots, and so forth can be defined with this system of exponentials.

For example, the square-root of (10^2) is clearly 10^1. That is to say: sqrt(11) == (10). The sqrt(10) == sqrt(10^1) == sqrt(10^4) (cycle property: just like how we can do this with "i") == 10^2 == 11.

sqrt(01) == sqrt(10^0) == 10^(0/2) == 10^0 == 01. sqrt(1) == 1 after all, in any and all well defined algebras.

-------------

Fractions also exist, as a unique inverse to any multiplication problem. For example: a / b == c, therefore, a = c * b^-1. You'll find that these numbers are consistent with their logarithms (!!!). Lets try out one example.

10 * 11 == 01

10 = 01 / 11 (aka: 10 is the multiplicative inverse of 11, or 10 == 11^-1. We also can see that 11 == 10^-1).

10^-1 == 10^2 (looping property of "i" as before), and as we have established earlier: 10^2 == 11.

---------

If it helps, maybe doing this in the complex field makes it more obvious:

(1/sqrt(2) * (1 + i)) ^ 1 == sqrt(2) * (1+i)

(1/sqrt(2) * (1 + i)) ^ 2 == 1/2 * (1 + 2i + i^2) == i

(1/sqrt(2) * (1 + i)) ^ 3 == i * 1/sqrt(2) * (1 + i) == 1/sqrt(2) * (i - 1)

(1/sqrt(2) * (1 + i)) ^ 4 == (1/sqrt(2) * (1 + i)) ^ 2 ^ 2 == i^2 == -1

And we can see that the cycle loops in the complex plane, albeit over cycle-8 instead of cycle-3, but the same concept.

A little Ot but, i questioned thinking, and humans, and some other stuff... with an result ^^

Solved: 'The Problem Is The Solution'

German Comic and hm it is written in a kind of dialect, so ask someone who might speak english (-:

Direct-Link to the ('Parted' Solution):

https://i.ibb.co/cFNMXGn/00422-NEIN-FINAL-Mail.png

https://i.ibb.co/TcCL856/00423-NEIN-FINAL-Mail.png

https://i.ibb.co/7YtGT0n/00424-NEIN-FINAL-Mail.png

regards, (-:

I wonder how does this affect deep learning calculations and how it's possible that all those multiplications and additions still work.

Reply on news.ycombinator.com