Re: improved norm behaviour

[Jaroslav Hajek]

On Mon, Aug 11, 2008 at 8:12 PM, John W. Eaton <address@hidden> wrote:
> On 11-Aug-2008, Jaroslav Hajek wrote:
>
> | gosh. Here it is.
> |
> |
> | On Mon, Aug 11, 2008 at 11:33 AM, Jaroslav Hajek <address@hidden> wrote:
> | > hello,
> | >
> | > please consider this changeset, that improves the builtin norm in
> | > several aspects.
> | > Summary of changes:
> | >
> | > 1. Everything is accessible from C++ (src/xnorm.h). __norm__ is no
> | > longer necessary (though not yet removed).
> | >
> | > 2. Scaled summation for vector p-norms to avoid overflows. Expressions
> | > such as `norm (1e20*ones (5,1), 20)' no longer produce an infinity.
> | >
> | > 3. norm (a, [p,] 'columns') or norm(a, [p,] 'rows') can compute column
> | > or row norms of a matrix. Currently, AFAIK there is no way in Octave
> | > (nor Matlab) to compute norms for a set of vectors both efficiently
> | > and robustly. Also, this operation is sometimes useful in linear
> | > algebra (scaling). Works efficiently on sparse matrices, too.
> | >
> | > 4. support for general matrix p-norms via Higham's hybrid algorithm.
> | > Seems to work well for real matrices, for complex matrices the
> | > accuracy can be poor. I intend to look into this eventually (the OSE
> | > step should probably be modified in complex case. The Higham's paper
> | > doesn't deal with complex cases, so it will take some investigation.)
> | >
> | > this patch won't work properly without my previous patch that fixes
> | > the segfaults produced by Array<T>::hermitian.
> | > Although 4. is partly incomplete, I think that the patch is useful
> | > enough to be included in this state. (I hope to contribute some tests
> | > eventually).
>
> What is the reason for eliminating the MArray<T>::norm functions?

Well, they lack the robustness robust (in case of general p-norm). I
wanted to alert any potential user that a new algorithm is available.

>
> Shouldn't a norm function always return a real-valued scalar?  If so,
> then wouldn't it be best to have most of the norm code in liboctave
> instead of liboctinterp, so that octave_value objects are not needed
> to compute the norm of Matrix, ColumnVector, or other numeric objects?
>

Good point. I guess I should move the non-octave_value to liboctave
(oct-norm.h + oct-norm.cc) I'd like to leave the xnorm.h and xnorm.cc
though, to keep the octave_value dispatchers out of data.cc, which is
already quite large. I'll update the patch and resubmit.

cheers,

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz