[math-fun] Floating point matrix product anomalies
From: Bill Gosper <billgosper@gmail.com> Subject: [math-fun] Floating point matrix product anomalies
--Wow, that is strange. Perhaps related is this: 1. Take N consecutive rightward steps of length 1+noise (where noise is a tiny-variance mean=0 random variable). You will reach N+-sqrt(N)*NoiseAmplitude. 2. Instead take a step, double it, add noise, and keep on doing that for lg(N) doublings. You will reach N+-N*Constant*NoiseAmplitude. 3. If in (1) the noise each step is negatively correlated to the noise the preceding step, then we do even better. -- Warren D. Smith http://RangeVoting.org <-- add your endorsement (by clicking "endorse" as 1st step)
* Warren D Smith <warren.wds@gmail.com> [Dec 07. 2015 16:16]:
From: Bill Gosper <billgosper@gmail.com> Subject: [math-fun] Floating point matrix product anomalies
--Wow, that is strange. Perhaps related is this:
1. Take N consecutive rightward steps of length 1+noise (where noise is a tiny-variance mean=0 random variable). You will reach N+-sqrt(N)*NoiseAmplitude.
2. Instead take a step, double it, add noise, and keep on doing that for lg(N) doublings. You will reach N+-N*Constant*NoiseAmplitude.
3. If in (1) the noise each step is negatively correlated to the noise the preceding step, then we do even better.
The example given by RWG is straight from the book of numeric cancellation (here 1 - cos(phi) for small values of phi). For a lack of a reference from the top of my head I point to Section 21.3.2 (p.417) of my fxtbook.
From the point of view of integration of differential equations: 1-step methods tend to not quite cut it.
About matrices: the one chosen can be diagonalized(*) as, say, A = U^{-1} diag( eigenvalues ) U so I'd compute A^n = U^-1 diag(eigenvalues)^n U. (*) it is anti-symmetric, so eigenvalues are purely imaginary and the eigenspace for every eigenvalue is full (geometric multiplicity == algebraic multiplicity for every eigenvalue). For U take the matrix whose columns are eigenvectors, no need to even normalize them, btw., as the term U^{-1} automagically takes care of that. If eigenvectors are a ONB (they can always be made to that) then U is orthonormal/unitary and U^{-1} = U^T (or U^{*} when actually complex). Best regards, jj
-- Warren D. Smith http://RangeVoting.org <-- add your endorsement (by clicking "endorse" as 1st step)
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com https://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
On 2015-12-07 07:34, Joerg Arndt wrote:
* Warren D Smith <warren.wds@gmail.com> [Dec 07. 2015 16:16]:
From: Bill Gosper <billgosper@gmail.com> Subject: [math-fun] Floating point matrix product anomalies
--Wow, that is strange. Perhaps related is this:
1. Take N consecutive rightward steps of length 1+noise (where noise is a tiny-variance mean=0 random variable). You will reach N+-sqrt(N)*NoiseAmplitude.
2. Instead take a step, double it, add noise, and keep on doing that for lg(N) doublings. You will reach N+-N*Constant*NoiseAmplitude.
3. If in (1) the noise each step is negatively correlated to the noise the preceding step, then we do even better.
The example given by RWG is straight from the book of numeric cancellation (here 1 - cos(phi) for small values of phi).
In what sense are pi/7 and 2 pi/7 small? --rwg
For a lack of a reference from the top of my head I point to Section 21.3.2 (p.417) of my fxtbook.
From the point of view of integration of differential equations: 1-step methods tend to not quite cut it.
About matrices: the one chosen can be diagonalized(*) as, say, A = U^{-1} diag( eigenvalues ) U so I'd compute A^n = U^-1 diag(eigenvalues)^n U. (*) it is anti-symmetric, so eigenvalues are purely imaginary and the eigenspace for every eigenvalue is full (geometric multiplicity == algebraic multiplicity for every eigenvalue). For U take the matrix whose columns are eigenvectors, no need to even normalize them, btw., as the term U^{-1} automagically takes care of that. If eigenvectors are a ONB (they can always be made to that) then U is orthonormal/unitary and U^{-1} = U^T (or U^{*} when actually complex).
Best regards, jj
-- Warren D. Smith
participants (3)
-
Joerg Arndt -
rwg -
Warren D Smith