bugGNU Octave - Bugs: bug #67610, precision of the octave norm...

 
 

bug #67610: precision of the octave norm function can be improved using the Kahan summation algorithm

Submitter:  Hendrik K <koerhen>
Submitted:  Thu 16 Oct 2025 09:50:50 AM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Patch Submitted Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Post a Comment

Add a New Comment Rich Markup
   

Discussion

Jump to the original submission

Fri 31 Oct 2025 07:31:25 PM UTC, comment #15: 

A similar discussion in MATLAB below.

https://blogs.mat ... ore-accurate-sum/

In short, since they split the sum among CPU cores there's already some built in protection.

Anonymous
Fri 24 Oct 2025 02:55:09 PM UTC, comment #14: 
Dmitri A. Sergatskov <dasergatskov>
Fri 24 Oct 2025 02:25:59 PM UTC, comment #13: 

The more I think about it, the more I convince that we should just use BLAS routines. We are re-inventing what is already done there.

WRT of Hammarling note -- perhaps you sent a wrong link, but I do not see how it is relevant to the code in oct-norm.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 24 Oct 2025 10:57:09 AM UTC, comment #12: 

And yes it is Hammarling and not Hammerling...



Hendrik K <koerhen>
Fri 24 Oct 2025 10:49:45 AM UTC, comment #11: 


@Hammerling: See
https://hammarlin ... lane-rots-IMA.pdf (1974)

and for an overview of the development in this space
https://perso.ens ... idesNormsPP22.pdf

So
i) Hammerling is timely BEFORE Blue (1978) AND
ii) Blue relies on three accumulators/scales and not just one like in oct-norm, so referencing Blue seems not to be correct when describing the algorithm used in oct-norm.

But please correct me if I am wrong.



Hendrik K <koerhen>
Fri 24 Oct 2025 10:48:48 AM UTC, comment #10: 

I like pairwise summation, but it is recursive from what I understand and the current SW architecture using accumulators with different algorithms/classes (inf norm, euclidean, 1-norm etc.) does ill support recursive algorithms or even the iterative, "shift-reduce" formulation for pairwise summation from wikipedia.

So not just the accumulator class for 1-norm would need changing, but also the calling (iterative) functions vector_norm, column_norms etc. when doing pairwise accumulation for norm-1. So it would be rather an intrusive change.

So what could be considered is implementing the faster BLAS summation approach (essentially building 3 buckets of running sums/accumulators: very small number, normal numbers, very large numbers) and sum the buckets together when the final result is requested. This would be likely as good as (and likely even better) than pairwise summation in term of precision and with much less overhead than Kahan and would fit very neatly in the current accumulator software structure.
The only drawbacks are
a) the algorithm is a more complex with the three running accumulator buckets than Kahan
b) the bucket definition/thresholds depends on the machine floating point implementation and data types (single/double). So a number of constants need to be determined (could be static variables) relying on e.g. https://en.cppref ... meric_limits.html, but I have not looked in detail.
 

So no real easy solution (either performance via Kahan or code complexity when solved differently as depicted above) for what I think is a more unlikely problem: Norm-one is the largest value by column, so for typical situations the sum is only for n values (and never for n*m values for an n*m matrix).
But yes - if someone would be willing to tackle this, it would improve Octave for norm-1 edge cases.


Side Note:
Bug report #61143 seems to be actually driven by two problems:
- spurious overflow (oct-norm uses scaling for this purpose; not required for 1-norm as spurious overflow it is impossible to occur in a simple sum) and
- precision issues (oct-norm will use Kahan for this purpose)
Both are solved using pairwise summation and mean building at the same time. #61143 do not have the issues octave would face in terms of existing SW architecture in oct-norm (iterative accumulative). But scaling or something else would be required on top for non-associative norms in oct-norm: Current scaling is an iterative algorithm, so it will pose a problem with any pairwise (recursive/vectorized) summation algorithm.

Hendrik K <koerhen>
Fri 24 Oct 2025 10:27:34 AM UTC, comment #9: 

I think "Hammerling" is a wrong attribution and it should be Blue's scaling (used in a BLAS routines NRM2 etc..) The reference would be:

Blue, J. L. (1978). "A Portable Fortran Program to Find the Euclidean Norm of a Vector". ACM Transactions on Mathematical Software, 4(1), 15-23.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 24 Oct 2025 09:39:44 AM UTC, comment #8: 

That may be a good idea. BTW, what is "Hammerling" scaling this code mentions? I suspect it could refer to Sven Hammarling (see https://hamma ... ing.com/sven/), but not avare of specific reference.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 24 Oct 2025 07:43:43 AM UTC, comment #7: 

I wonder if we should think about pairwise summation as an intermediate compromise between naive summation and full Kahan algorithm.

Bug report #61143 has a long and worthy discussion of the topic.

Also see Wikipedia entry on subject at https://en.wikipe ... airwise_summation.



Rik <rik5>
Group administrator
Sun 19 Oct 2025 01:56:39 AM UTC, comment #6: 
Anonymous
Thu 16 Oct 2025 05:15:27 PM UTC, comment #5: 

Apparently on ARM64 both Netlib and OpenBlas NRM2 is imprecise.
I guess this is another vote for using kahan summation (for now).
I still think that `norm-1` should be modified to use it as well.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 16 Oct 2025 01:05:22 PM UTC, comment #4: 

For me this argument works the other way -- what is the point of having some function if it is just basic sum? When I call norm I expect it to calcalates norm and check some edge cases...

Also I suspect matlab (not sure about Julia) does not calculate sum as a straight incremented accumulator, we had this discussion many times already, it definitely does some kind of compensation (though probably not Kahan).

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 16 Oct 2025 11:46:44 AM UTC, comment #3: 

Typo removed, white spaces removed and indent of if then else with curly brackets changed and update patch attached.



@norm-1:

I did not do benchmarks, but
a)  we would replace just one (basic) operation/line:
 m_sum += std::abs (val);

with several basic operations/lines when using Kahan. So this will result in a high relative overhead.

b) On top norm 1 just sums the individual (absolute) values by each column without any modification. In contrast to this is norm2: pow2  typically explodes the range of values making very small and very large values leading to very likely summation potential numerical issues.
Stated differently: we are doing kind of a simple sum(x). Do we use Kahan for the normal octave sum function as well ? (Small side note: Neither Matlab nor Julia use Kahan for sum - I guess for good reasons...)

(file #57724)

Hendrik K <koerhen>
Thu 16 Oct 2025 10:39:04 AM UTC, comment #2: 

There are some style issues:
"the algorithms uses" ->  "the algorithm uses"

Else-if loop braces should be

else if (t != 0)
  {
    R y = std::pow (t/m_scl, m_p) - m_comp;
    ...
  }

Some trailing spaces here and there.

---------

I am also not convinced that Kahan should not be used for norm-1.
Do you have actual benchamarks?

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 16 Oct 2025 09:59:00 AM UTC, comment #1: 

Patch is attached.

(file #57723)

Hendrik K <koerhen>
Thu 16 Oct 2025 09:50:50 AM UTC, original submission:  

The octave norm function uses an accumulation algorithm which can numerically be improved:
 
Currently it uses scaling (simple Hammerling approach) to avoid spurious overflows and underflows and then simply adds the individual (scaled) terms to the running sum.
 
Using the more sophisticated Kahan algorithm for summing up, the numerical precision can significantly be improved (in exchange of some performance impacts due to the required higher number of numerical operations).

All accumulators now use the Kahan algorithm except the norm-1 accumulator (the potential numerical accumulation error is materially smaller here and the relative Kahan performance overhead is materially higher: norm-1 is simply the the largest column sum of the absolute values of A, so no scaling and pow2/std::pow required)

See https://octave.di ... in-octave/7034/16 for a more detailed discussion.

A patch is attached and make check is passed.


Hendrik K <koerhen>

 

Attached Files

Attached Files
file #57724:  67610_Update.patch added by koerhen (3.0KiB - text/x-patch)
file #57723:  67610.patch added by koerhen (3.1KiB - text/x-patch)

(Note: upload size limit is set to 16MiB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

 

Dependencies

This item does not depend on any other items.

No items depend on this one.

 

Mail Notification Carbon-Copy List

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by koerhen (Submitted the item)
  •  

    Votes

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    History

    Follow 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2025-10-24 rik5 StatusNone Patch Submitted
    2025-10-16 koerhen Attached File- Added 67610_Update.patch, #57724
    2025-10-16 koerhen Attached File- Added 67610.patch, #57723

    Back to the top

    Powered by Savane 3.16.
    Corresponding source code