20

Intel support for IEEE float16 storage format Intel supports IEEE half as a storage type in processors since Ivy Bridge (2013). Storage type means you can get a memory/cache capacity/bandwidth advantage but the compute is done with single precision after converting to and from the IEEE half precision format. https://software.intel.com/content/www/us/en/...


7

In my opinion, not very uniformly. Low precision arithmetic seems to have gained some traction in machine learning, but there's varying definitions for what people mean by low precision. There's the IEEE-754 half (10 bit mantissa, 5 bit exponent, 1 bit sign) but also bfloat16 (7 bit mantissa, 8 bit exponent, 1 bit sign) which favors dynamic range over ...


6

The accepted answer provides an overview. I'll add a few more details about support in NVIDIA processors. The support I'm describing here is 16 bit, IEEE 754 compliant, floating point arithmetic support, including add, multiply, multiply-add, and conversions to/from other formats. Maxwell (circa 2015) The earliest IEEE 754 FP16 ("binary16" or "half ...


5

Orthogonal matrices are about as well-conditioned as you can get, but numerical errors still occur. One common error is loss of orthogonality. A fix for this could be to re-orthogonalize your columns after some number of multiplications. You can do this by just taking the QR decomposition of your matrix after some number of products and taking the orthogonal ...


5

There is no need for numerical computation here. First, $T(q)$ is a well-known function, the logarithmic integral. Repeated integration by parts gives an asymptotic expansion $$\mathrm{Li}(q) = \frac{q}{\log q}\sum_{k=0}^{K-1} \frac{k!}{\log^k q} + O\left(\frac{q}{\log^{K+1}q}\right).$$ There's also a fairly rapidly convergent representation due to ...


4

You can find out if your hardware supports half-precision via: $ lscpu | grep Flags | grep f16c Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf ...


4

You write $S(q)$ and $T(q)$ as integrals, but it is easier to think of them as solutions of ODEs: $$ S'(q) = \sin^2\left(\frac{π\Gamma(q)}{2q}\right) $$ with initial conditions $$ S(2) = 0, $$ and similarly for $T(q)$. You can then use any of the common ODE integrators in matlab, mathematica, maple, ..., to solve and plot the solutions so that you can ...


4

Very interesting question! LAPACK-inspired adaptive strategy This reminds me of a bug that was found in a LAPACK routine (rank-revealing QR) related to 'downdating' norms: essentially, you are given the norm of a vector v, and you want to compute at each iteration in O(1) the norm of the same vector after chopping off its initial entry: v[1:], v[2:], ... (...


3

You don't need the inverse even with the goal of finding $K^{-1} h h^{T} K^{-1} - K^{-1}$. If you are interested to have this expression, I would explain how you can convert it to a matrix equation and then solve it more efficiently: Let's define the $X$ as: $$X = K^{-1} h h^{T} K^{-1} - K^{-1}$$ Your objective is to calculate $X$ in this equation by ...


3

After thinking about this some more, I can answer this one myself! I don't think the complex plane makes the log-sum-exp trick appreciably different, at least in Cartesian coordinates. In particular, if $z=u+iv$ then $e^z=e^{u+iv}=e^u (\cos v + i\sin v).$ Notice the $v$ part has magnitude 1 by construction, so overflow or underflow is principally caused ...


3

In this answer, I recommended using mpmath Python library for arbitrary precision. Notice, that since matrices in mpmath are implemented as dictionaries: Only non-zero values are stored, so it is cheap to represent sparse matrices. thus, this particular library seems like a good fit for your purpose of debugging. I find this library (and approach in ...


3

The only chance you stand to deal with this problem from a numerical perspective is oscillatory integration methods. Filon/Levin-type methods can sometimes handle problems like this, particularly when they are of $\sin$ or $\cos$ type, though the $\Gamma$ function and its run-away growth may be prohibitive for the $q$ you are hoping for. In any case, I was ...


3

This particular question is related to the fields theory of computation and computational mathematics(or scientific computing), but some of aspects of it fall under computer science. Nick Higham's "Accuracy and Stability of Numerical Algorithms" is a must-read. Though while the first few chapters are fairly general, the most of the book focuses on ...


2

If you first convert to non-normalized form and then normalize, you will of course get the zeros at the end simply because the normalization can not know what else to put in there. But that's not how one converts decimal numbers into floating point. One does the normalization first, of course, and only then computes the digits past the (binary) point.


2

I think the primary issue you're seeing is related to how you're checking for convergence. The chosen value for your tolerance is likely too tight. The tolerance should at a minimum allow 1 double precision (DP) epsilon up or down of floating point rounding error. However, 1 DP epsilon is relative to the magnitude of the floating point value. In this ...


1

Based on Federico's suggestions and ideas, more straight forward formulation of extracting $K^{-1}hh^{T}K^{-1}-K^{-1}$ would be: $$X = K^{-1}hh^{T}K^{-1}-K^{-1}$$ $$KXK = hh^{T}-K$$ $$Z = XK$$ Solve for $Z$: $$KZ = hh^{T} - K$$ and then find $X^{T}$ from: $$K^{T} X^{T} = Z^{T}$$ and finally $X$: $$X = (X^{T})^{T}$$ Let's define the error between ...


1

A simpler scheme In addition to my other answer which is more about an extension of the Fornberg method, I'll address here the question for more simple alternatives. For this I sketch an alternative scheme which produces the derivative coefficients the of Lagrangian interpolation more directly. Its implementation requires only a few lines of code, works ...


Only top voted, non community-wiki answers of a minimum length are eligible