oct-norm.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 2008-2012 VZLU Prague, a.s.
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 // author: Jaroslav Hajek <highegg@gmail.com>
00024 
00025 #ifdef HAVE_CONFIG_H
00026 #include <config.h>
00027 #endif
00028 
00029 #include <cassert>
00030 #include <cfloat>
00031 #include <cmath>
00032 
00033 #include <iostream>
00034 #include <vector>
00035 
00036 #include "oct-cmplx.h"
00037 #include "lo-error.h"
00038 #include "lo-ieee.h"
00039 #include "Array.h"
00040 #include "Array-util.h"
00041 #include "CMatrix.h"
00042 #include "dMatrix.h"
00043 #include "fCMatrix.h"
00044 #include "fMatrix.h"
00045 #include "CColVector.h"
00046 #include "dColVector.h"
00047 #include "CRowVector.h"
00048 #include "dRowVector.h"
00049 #include "fCColVector.h"
00050 #include "fColVector.h"
00051 #include "fCRowVector.h"
00052 #include "fRowVector.h"
00053 #include "CSparse.h"
00054 #include "dSparse.h"
00055 #include "dbleSVD.h"
00056 #include "CmplxSVD.h"
00057 #include "floatSVD.h"
00058 #include "fCmplxSVD.h"
00059 
00060 // Theory: norm accumulator is an object that has an accum method able
00061 // to handle both real and complex element, and a cast operator
00062 // returning the intermediate norm. Reference: Higham, N. "Estimating
00063 // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992.
00064 
00065 // norm accumulator for the p-norm
00066 template <class R>
00067 class norm_accumulator_p
00068 {
00069   R p,scl,sum;
00070 public:
00071   norm_accumulator_p () {} // we need this one for Array
00072   norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {}
00073 
00074   template<class U>
00075   void accum (U val)
00076     {
00077       octave_quit ();
00078       R t = std::abs (val);
00079       if (scl == t) // we need this to handle Infs properly
00080         sum += 1;
00081       else if (scl < t)
00082         {
00083           sum *= std::pow (scl/t, p);
00084           sum += 1;
00085           scl = t;
00086         }
00087       else if (t != 0)
00088         sum += std::pow (t/scl, p);
00089     }
00090   operator R () { return scl * std::pow (sum, 1/p); }
00091 };
00092 
00093 // norm accumulator for the minus p-pseudonorm
00094 template <class R>
00095 class norm_accumulator_mp
00096 {
00097   R p,scl,sum;
00098 public:
00099   norm_accumulator_mp () {} // we need this one for Array
00100   norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {}
00101 
00102   template<class U>
00103   void accum (U val)
00104     {
00105       octave_quit ();
00106       R t = 1 / std::abs (val);
00107       if (scl == t)
00108         sum += 1;
00109       else if (scl < t)
00110         {
00111           sum *= std::pow (scl/t, p);
00112           sum += 1;
00113           scl = t;
00114         }
00115       else if (t != 0)
00116         sum += std::pow (t/scl, p);
00117     }
00118   operator R () { return scl * std::pow (sum, -1/p); }
00119 };
00120 
00121 // norm accumulator for the 2-norm (euclidean)
00122 template <class R>
00123 class norm_accumulator_2
00124 {
00125   R scl,sum;
00126   static R pow2 (R x) { return x*x; }
00127 public:
00128   norm_accumulator_2 () : scl(0), sum(1) {}
00129 
00130   void accum (R val)
00131     {
00132       R t = std::abs (val);
00133       if (scl == t)
00134         sum += 1;
00135       else if (scl < t)
00136         {
00137           sum *= pow2 (scl/t);
00138           sum += 1;
00139           scl = t;
00140         }
00141       else if (t != 0)
00142         sum += pow2 (t/scl);
00143     }
00144 
00145   void accum (std::complex<R> val)
00146     {
00147       accum (val.real ());
00148       accum (val.imag ());
00149     }
00150 
00151   operator R () { return scl * std::sqrt (sum); }
00152 };
00153 
00154 // norm accumulator for the 1-norm (city metric)
00155 template <class R>
00156 class norm_accumulator_1
00157 {
00158   R sum;
00159 public:
00160   norm_accumulator_1 () : sum (0) {}
00161   template<class U>
00162   void accum (U val)
00163     {
00164       sum += std::abs (val);
00165     }
00166   operator R () { return sum; }
00167 };
00168 
00169 // norm accumulator for the inf-norm (max metric)
00170 template <class R>
00171 class norm_accumulator_inf
00172 {
00173   R max;
00174 public:
00175   norm_accumulator_inf () : max (0) {}
00176   template<class U>
00177   void accum (U val)
00178     {
00179       max = std::max (max, std::abs (val));
00180     }
00181   operator R () { return max; }
00182 };
00183 
00184 // norm accumulator for the -inf pseudonorm (min abs value)
00185 template <class R>
00186 class norm_accumulator_minf
00187 {
00188   R min;
00189 public:
00190   norm_accumulator_minf () : min (octave_Inf) {}
00191   template<class U>
00192   void accum (U val)
00193     {
00194       min = std::min (min, std::abs (val));
00195     }
00196   operator R () { return min; }
00197 };
00198 
00199 // norm accumulator for the 0-pseudonorm (hamming distance)
00200 template <class R>
00201 class norm_accumulator_0
00202 {
00203   unsigned int num;
00204 public:
00205   norm_accumulator_0 () : num (0) {}
00206   template<class U>
00207   void accum (U val)
00208     {
00209       if (val != static_cast<U> (0)) ++num;
00210     }
00211   operator R () { return num; }
00212 };
00213 
00214 
00215 // OK, we're armed :) Now let's go for the fun
00216 
00217 template <class T, class R, class ACC>
00218 inline void vector_norm (const Array<T>& v, R& res, ACC acc)
00219 {
00220   for (octave_idx_type i = 0; i < v.numel (); i++)
00221     acc.accum (v(i));
00222 
00223   res = acc;
00224 }
00225 
00226 // dense versions
00227 template <class T, class R, class ACC>
00228 void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
00229 {
00230   res = MArray<R> (dim_vector (1, m.columns ()));
00231   for (octave_idx_type j = 0; j < m.columns (); j++)
00232     {
00233       ACC accj = acc;
00234       for (octave_idx_type i = 0; i < m.rows (); i++)
00235         accj.accum (m(i, j));
00236 
00237       res.xelem (j) = accj;
00238     }
00239 }
00240 
00241 template <class T, class R, class ACC>
00242 void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
00243 {
00244   res = MArray<R> (dim_vector (m.rows (), 1));
00245   std::vector<ACC> acci (m.rows (), acc);
00246   for (octave_idx_type j = 0; j < m.columns (); j++)
00247     {
00248       for (octave_idx_type i = 0; i < m.rows (); i++)
00249         acci[i].accum (m(i, j));
00250     }
00251 
00252   for (octave_idx_type i = 0; i < m.rows (); i++)
00253     res.xelem (i) = acci[i];
00254 }
00255 
00256 // sparse versions
00257 template <class T, class R, class ACC>
00258 void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
00259 {
00260   res = MArray<R> (dim_vector (1, m.columns ()));
00261   for (octave_idx_type j = 0; j < m.columns (); j++)
00262     {
00263       ACC accj = acc;
00264       for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
00265         accj.accum (m.data (k));
00266 
00267       res.xelem (j) = accj;
00268     }
00269 }
00270 
00271 template <class T, class R, class ACC>
00272 void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
00273 {
00274   res = MArray<R> (dim_vector (m.rows (), 1));
00275   std::vector<ACC> acci (m.rows (), acc);
00276   for (octave_idx_type j = 0; j < m.columns (); j++)
00277     {
00278       for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
00279         acci[m.ridx (k)].accum (m.data (k));
00280     }
00281 
00282   for (octave_idx_type i = 0; i < m.rows (); i++)
00283     res.xelem (i) = acci[i];
00284 }
00285 
00286 // now the dispatchers
00287 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \
00288 template <class T, class R> \
00289 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \
00290 { \
00291   RES_TYPE res; \
00292   if (p == 2) \
00293     FUNC_NAME (v, res, norm_accumulator_2<R> ()); \
00294   else if (p == 1) \
00295     FUNC_NAME (v, res, norm_accumulator_1<R> ()); \
00296   else if (lo_ieee_isinf (p)) \
00297     { \
00298       if (p > 0) \
00299         FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \
00300       else \
00301         FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \
00302     } \
00303   else if (p == 0) \
00304     FUNC_NAME (v, res, norm_accumulator_0<R> ()); \
00305   else if (p > 0) \
00306     FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \
00307   else \
00308     FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \
00309   return res; \
00310 }
00311 
00312 DEFINE_DISPATCHER (vector_norm, MArray<T>, R)
00313 DEFINE_DISPATCHER (column_norms, MArray<T>, MArray<R>)
00314 DEFINE_DISPATCHER (row_norms, MArray<T>, MArray<R>)
00315 DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>)
00316 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>)
00317 
00318 // The approximate subproblem in Higham's method. Find lambda and mu such that
00319 // norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized.
00320 // Real version. As in Higham's paper.
00321 template <class ColVectorT, class R>
00322 static void
00323 higham_subp (const ColVectorT& y, const ColVectorT& col,
00324              octave_idx_type nsamp, R p, R& lambda, R& mu)
00325 {
00326   R nrm = 0;
00327   for (octave_idx_type i = 0; i < nsamp; i++)
00328     {
00329       octave_quit ();
00330       R fi = i*M_PI/nsamp, lambda1 = cos (fi), mu1 = sin (fi);
00331       R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
00332                          std::pow (std::abs (mu1), p), 1/p);
00333       lambda1 /= lmnr; mu1 /= lmnr;
00334       R nrm1 = vector_norm (lambda1 * y + mu1 * col, p);
00335       if (nrm1 > nrm)
00336         {
00337           lambda = lambda1;
00338           mu = mu1;
00339           nrm = nrm1;
00340         }
00341     }
00342 }
00343 
00344 // Complex version. Higham's paper does not deal with complex case, so we use a simple
00345 // extension. First, guess the magnitudes as in real version, then try to rotate lambda
00346 // to improve further.
00347 template <class ColVectorT, class R>
00348 static void
00349 higham_subp (const ColVectorT& y, const ColVectorT& col,
00350              octave_idx_type nsamp, R p,
00351              std::complex<R>& lambda, std::complex<R>& mu)
00352 {
00353   typedef std::complex<R> CR;
00354   R nrm = 0;
00355   lambda = 1.0;
00356   CR lamcu = lambda / std::abs (lambda);
00357   // Probe magnitudes
00358   for (octave_idx_type i = 0; i < nsamp; i++)
00359     {
00360       octave_quit ();
00361       R fi = i*M_PI/nsamp, lambda1 = cos (fi), mu1 = sin (fi);
00362       R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
00363                          std::pow (std::abs (mu1), p), 1/p);
00364       lambda1 /= lmnr; mu1 /= lmnr;
00365       R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p);
00366       if (nrm1 > nrm)
00367         {
00368           lambda = lambda1 * lamcu;
00369           mu = mu1;
00370           nrm = nrm1;
00371         }
00372     }
00373   R lama = std::abs (lambda);
00374   // Probe orientation
00375   for (octave_idx_type i = 0; i < nsamp; i++)
00376     {
00377       octave_quit ();
00378       R fi = i*M_PI/nsamp;
00379       lamcu = CR (cos (fi), sin (fi));
00380       R nrm1 = vector_norm (lama * lamcu * y + mu * col, p);
00381       if (nrm1 > nrm)
00382         {
00383           lambda = lama * lamcu;
00384           nrm = nrm1;
00385         }
00386     }
00387 }
00388 
00389 // the p-dual element (should work for both real and complex)
00390 template <class T, class R>
00391 inline T elem_dual_p (T x, R p)
00392 {
00393   return signum (x) * std::pow (std::abs (x), p-1);
00394 }
00395 
00396 // the VectorT is used for vectors, but actually it has to be
00397 // a Matrix type to allow all the operations. For instance SparseMatrix
00398 // does not support multiplication with column/row vectors.
00399 // the dual vector
00400 template <class VectorT, class R>
00401 VectorT dual_p (const VectorT& x, R p, R q)
00402 {
00403   VectorT res (x.dims ());
00404   for (octave_idx_type i = 0; i < x.numel (); i++)
00405     res.xelem (i) = elem_dual_p (x(i), p);
00406   return res / vector_norm (res, q);
00407 }
00408 
00409 // Higham's hybrid method
00410 template <class MatrixT, class VectorT, class R>
00411 R higham (const MatrixT& m, R p, R tol, int maxiter,
00412           VectorT& x)
00413 {
00414   x.resize (m.columns (), 1);
00415   // the OSE part
00416   VectorT y(m.rows (), 1, 0), z(m.rows (), 1);
00417   typedef typename VectorT::element_type RR;
00418   RR lambda = 0, mu = 0;
00419   for (octave_idx_type k = 0; k < m.columns (); k++)
00420     {
00421       octave_quit ();
00422       VectorT col (m.column (k));
00423       if (k > 0)
00424         higham_subp (y, col, 4*k, p, lambda, mu);
00425       for (octave_idx_type i = 0; i < k; i++)
00426         x(i) *= lambda;
00427       x(k) = mu;
00428       y = lambda * y + mu * col;
00429     }
00430 
00431   // the PM part
00432   x = x / vector_norm (x, p);
00433   R q = p/(p-1);
00434 
00435   R gamma = 0, gamma1;
00436   int iter = 0;
00437   while (iter < maxiter)
00438     {
00439       octave_quit ();
00440       y = m*x;
00441       gamma1 = gamma;
00442       gamma = vector_norm (y, p);
00443       z = dual_p (y, p, q);
00444       z = z.hermitian ();
00445       z = z * m;
00446 
00447       if (iter > 0 && (vector_norm (z, q) <= gamma
00448                        || (gamma - gamma1) <= tol*gamma))
00449         break;
00450 
00451       z = z.hermitian ();
00452       x = dual_p (z, q, p);
00453       iter ++;
00454     }
00455 
00456   return gamma;
00457 }
00458 
00459 // derive column vector and SVD types
00460 
00461 static const char *p_less1_gripe = "xnorm: p must be at least 1";
00462 
00463 // Static constant to control the maximum number of iterations. 100 seems to be a good value.
00464 // Eventually, we can provide a means to change this constant from Octave.
00465 static int max_norm_iter = 100;
00466 
00467 // version with SVD for dense matrices
00468 template <class MatrixT, class VectorT, class SVDT, class R>
00469 R matrix_norm (const MatrixT& m, R p, VectorT, SVDT)
00470 {
00471   R res = 0;
00472   if (p == 2)
00473     {
00474       SVDT svd (m, SVD::sigma_only);
00475       res = svd.singular_values () (0,0);
00476     }
00477   else if (p == 1)
00478     res = xcolnorms (m, 1).max ();
00479   else if (lo_ieee_isinf (p))
00480     res = xrownorms (m, 1).max ();
00481   else if (p > 1)
00482     {
00483       VectorT x;
00484       const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
00485       res = higham (m, p, sqrteps, max_norm_iter, x);
00486     }
00487   else
00488     (*current_liboctave_error_handler) (p_less1_gripe);
00489 
00490   return res;
00491 }
00492 
00493 // SVD-free version for sparse matrices
00494 template <class MatrixT, class VectorT, class R>
00495 R matrix_norm (const MatrixT& m, R p, VectorT)
00496 {
00497   R res = 0;
00498   if (p == 1)
00499     res = xcolnorms (m, 1).max ();
00500   else if (lo_ieee_isinf (p))
00501     res = xrownorms (m, 1).max ();
00502   else if (p > 1)
00503     {
00504       VectorT x;
00505       const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
00506       res = higham (m, p, sqrteps, max_norm_iter, x);
00507     }
00508   else
00509     (*current_liboctave_error_handler) (p_less1_gripe);
00510 
00511   return res;
00512 }
00513 
00514 // and finally, here's what we've promised in the header file
00515 
00516 #define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \
00517   OCTAVE_API RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \
00518   { return vector_norm (x, p); } \
00519   OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \
00520   { return vector_norm (x, p); } \
00521   OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \
00522   { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \
00523   OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \
00524   { return vector_norm (x, static_cast<RTYPE> (2)); }
00525 
00526 DEFINE_XNORM_FUNCS(, double)
00527 DEFINE_XNORM_FUNCS(Complex, double)
00528 DEFINE_XNORM_FUNCS(Float, float)
00529 DEFINE_XNORM_FUNCS(FloatComplex, float)
00530 
00531 // this is needed to avoid copying the sparse matrix for xfrobnorm
00532 template <class T, class R>
00533 inline void array_norm_2 (const T* v, octave_idx_type n, R& res)
00534 {
00535   norm_accumulator_2<R> acc;
00536   for (octave_idx_type i = 0; i < n; i++)
00537     acc.accum (v[i]);
00538 
00539   res = acc;
00540 }
00541 
00542 #define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \
00543   OCTAVE_API RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \
00544   { return matrix_norm (x, p, PREFIX##Matrix ()); } \
00545   OCTAVE_API RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \
00546   { \
00547     RTYPE res; \
00548     array_norm_2 (x.data (), x.nnz (), res); \
00549     return res; \
00550   }
00551 
00552 DEFINE_XNORM_SPARSE_FUNCS(, double)
00553 DEFINE_XNORM_SPARSE_FUNCS(Complex, double)
00554 
00555 #define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \
00556   extern OCTAVE_API RPREFIX##RowVector xcolnorms (const PREFIX##Matrix& m, RTYPE p) \
00557   { return column_norms (m, p); } \
00558   extern OCTAVE_API RPREFIX##ColumnVector xrownorms (const PREFIX##Matrix& m, RTYPE p) \
00559   { return row_norms (m, p); } \
00560 
00561 DEFINE_COLROW_NORM_FUNCS(, , double)
00562 DEFINE_COLROW_NORM_FUNCS(Complex, , double)
00563 DEFINE_COLROW_NORM_FUNCS(Float, Float, float)
00564 DEFINE_COLROW_NORM_FUNCS(FloatComplex, Float, float)
00565 
00566 DEFINE_COLROW_NORM_FUNCS(Sparse, , double)
00567 DEFINE_COLROW_NORM_FUNCS(SparseComplex, , double)
00568 
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines