Boost logo

Glas :

Re: [glas] norms and inner products for vectors of matrices

From: Peter Gottschling (pgottsch_at_[hidden])
Date: 2006-01-13 10:41:58


On 13.01.2006, at 05:02, Karl Meerbergen wrote:

> Peter Gottschling wrote:
>
>> Concerning the norms, I would like to limit the result type to real
>> values as I have seen it in all mathematical definitions. Does
>> somebody sees a reason to vectors or matrices as results?
>>
>> The definitions of vector<vector> norms are straight forward:
>> - norm_1 := sum (norm_1(x[i])
>> - norm_2 := sqrt(sum(norm_2(x[i])^2)) // there might be
>> more
>> efficient ways to compute this
>> - norm_inf := max(norm_inf(x[i])
>>
>> Best,
>> Peter
>>
>>
>>
>>
>
> I am afraid that we cannot define norm_1, etc in this way. (I too made
> this mistake in a previous mail) But we can define functions
> generalized_norm_1 := sum (generalized_norm_1(x[i])
> generalized_norm_2 :=sqrt( sum (generalized_norm_2(x[i])^2 )
> generalized_norm_inf := max (generalized_norm_inf(x[i])
>
> norm_1, norm_2, norm_inf are well defined in linear algebra and these
> definitions should be respected:
> norm_1 := sum( abs(x[i]) )
> etc.
> These should only be used for vector and matrix objects with 'scalar'
> value_type's.
>
>
> Karl
>
Karl, I agree that it is not a good idea to call always the recursive
definition.

(Even if I think it is possible by defining norm_1 of scalars as abs
etc. Anyway, it's ugly.)

The way I would like to define it is to have one function and dispatch
first between scalar value types and nested value types. For instance

template <typename Vector>
double norm_1(const Vector& v, type::scalar)
{
     return sum(abs(x[i])); // well this is only pseudo-code
}

template <typename Vector>
double norm_1(const Vector& v, type::nested)
{
     return sum(norm_1(x[i]));
}

template <typename Vector>
double norm_1(const Vector& v)
{
     nesting_traits<typename Vector::value_type>::type nesting;
     return norm_1(v, nesting);
}

With this code we have the classical defintion if the value_type is
scalar and the recursive definition is the value_type itself is a
vector. Note that this also handles higher levels of nesting than 2.

What remains is the question how we define norm_1(vector<matrix<T> >) ?
  I'm not convinced that my definition is perfect for this. OTOH, I
wonder if there is a reasonable norm definition for vector<matrix> ????

May be we should discuss first what we consider the norm of a matrix.
Matrices are typically used as linear operators and operator norms can
be defined in terms of vector norms norm_xx(op)=sup( norm_xx( op(v) ) /
norm_xx (v) ). I vote for defining matrix norms as operator norms.
OTOH, matrices model the concept of VectorSpace and it is
mathematically perfectly correct to apply the vector norm directly. To
do this I suggest to use views for it, e.g.

vector_view<matrix<...> > vv(m);
alpha= norm_xx(vv);

One complication I see is that nested matrices has to be considered as
nested vectors, in other words the view must be applied recursively.

Coming back to the question what is norm_xx( vector<matrix> ) it would
be a mixture of vector and operator norms, which doesn't look
reasonable to me. Probably we see it clearer if we know what usage
vector<matrix> have (if there is any).

Best,
Peter
------------
Peter Gottschling
Research Associate
Open Systems Laboratory
Indiana University
301i Lindley Hall
Bloomington, IN 47405
Tel.: +1 812 855-8898 Fax: +1 812 856 0853
http://www.osl.iu.edu/~pgottsch