GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
oct-norm.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2008-2015 VZLU Prague, a.s.
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 // author: Jaroslav Hajek <highegg@gmail.com>
24 
25 #ifdef HAVE_CONFIG_H
26 #include <config.h>
27 #endif
28 
29 #include <cassert>
30 #include <cfloat>
31 #include <cmath>
32 
33 #include <iostream>
34 #include <vector>
35 
36 #include "oct-cmplx.h"
37 #include "lo-error.h"
38 #include "lo-ieee.h"
39 #include "mx-cm-s.h"
40 #include "mx-s-cm.h"
41 #include "mx-fcm-fs.h"
42 #include "mx-fs-fcm.h"
43 #include "Array.h"
44 #include "Array-util.h"
45 #include "CMatrix.h"
46 #include "dMatrix.h"
47 #include "fCMatrix.h"
48 #include "fMatrix.h"
49 #include "CColVector.h"
50 #include "dColVector.h"
51 #include "CRowVector.h"
52 #include "dRowVector.h"
53 #include "fCColVector.h"
54 #include "fColVector.h"
55 #include "fCRowVector.h"
56 #include "fRowVector.h"
57 #include "CSparse.h"
58 #include "dSparse.h"
59 #include "dbleSVD.h"
60 #include "CmplxSVD.h"
61 #include "floatSVD.h"
62 #include "fCmplxSVD.h"
63 
64 // Theory: norm accumulator is an object that has an accum method able
65 // to handle both real and complex element, and a cast operator
66 // returning the intermediate norm. Reference: Higham, N. "Estimating
67 // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992.
68 
69 // norm accumulator for the p-norm
70 template <class R>
72 {
73  R p,scl,sum;
74 public:
75  norm_accumulator_p () {} // we need this one for Array
76  norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {}
77 
78  template<class U>
79  void accum (U val)
80  {
81  octave_quit ();
82  R t = std::abs (val);
83  if (scl == t) // we need this to handle Infs properly
84  sum += 1;
85  else if (scl < t)
86  {
87  sum *= std::pow (scl/t, p);
88  sum += 1;
89  scl = t;
90  }
91  else if (t != 0)
92  sum += std::pow (t/scl, p);
93  }
94  operator R () { return scl * std::pow (sum, 1/p); }
95 };
96 
97 // norm accumulator for the minus p-pseudonorm
98 template <class R>
100 {
101  R p,scl,sum;
102 public:
103  norm_accumulator_mp () {} // we need this one for Array
104  norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {}
105 
106  template<class U>
107  void accum (U val)
108  {
109  octave_quit ();
110  R t = 1 / std::abs (val);
111  if (scl == t)
112  sum += 1;
113  else if (scl < t)
114  {
115  sum *= std::pow (scl/t, p);
116  sum += 1;
117  scl = t;
118  }
119  else if (t != 0)
120  sum += std::pow (t/scl, p);
121  }
122  operator R () { return scl * std::pow (sum, -1/p); }
123 };
124 
125 // norm accumulator for the 2-norm (euclidean)
126 template <class R>
128 {
129  R scl,sum;
130  static R pow2 (R x) { return x*x; }
131 public:
132  norm_accumulator_2 () : scl(0), sum(1) {}
133 
134  void accum (R val)
135  {
136  R t = std::abs (val);
137  if (scl == t)
138  sum += 1;
139  else if (scl < t)
140  {
141  sum *= pow2 (scl/t);
142  sum += 1;
143  scl = t;
144  }
145  else if (t != 0)
146  sum += pow2 (t/scl);
147  }
148 
149  void accum (std::complex<R> val)
150  {
151  accum (val.real ());
152  accum (val.imag ());
153  }
154 
155  operator R () { return scl * std::sqrt (sum); }
156 };
157 
158 // norm accumulator for the 1-norm (city metric)
159 template <class R>
161 {
162  R sum;
163 public:
164  norm_accumulator_1 () : sum (0) {}
165  template<class U>
166  void accum (U val)
167  {
168  sum += std::abs (val);
169  }
170  operator R () { return sum; }
171 };
172 
173 // norm accumulator for the inf-norm (max metric)
174 template <class R>
176 {
177  R max;
178 public:
179  norm_accumulator_inf () : max (0) {}
180  template<class U>
181  void accum (U val)
182  {
183  max = std::max (max, std::abs (val));
184  }
185  operator R () { return max; }
186 };
187 
188 // norm accumulator for the -inf pseudonorm (min abs value)
189 template <class R>
191 {
192  R min;
193 public:
195  template<class U>
196  void accum (U val)
197  {
198  min = std::min (min, std::abs (val));
199  }
200  operator R () { return min; }
201 };
202 
203 // norm accumulator for the 0-pseudonorm (hamming distance)
204 template <class R>
206 {
207  unsigned int num;
208 public:
209  norm_accumulator_0 () : num (0) {}
210  template<class U>
211  void accum (U val)
212  {
213  if (val != static_cast<U> (0)) ++num;
214  }
215  operator R () { return num; }
216 };
217 
218 
219 // OK, we're armed :) Now let's go for the fun
220 
221 template <class T, class R, class ACC>
222 inline void vector_norm (const Array<T>& v, R& res, ACC acc)
223 {
224  for (octave_idx_type i = 0; i < v.numel (); i++)
225  acc.accum (v(i));
226 
227  res = acc;
228 }
229 
230 // dense versions
231 template <class T, class R, class ACC>
232 void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
233 {
234  res = MArray<R> (dim_vector (1, m.columns ()));
235  for (octave_idx_type j = 0; j < m.columns (); j++)
236  {
237  ACC accj = acc;
238  for (octave_idx_type i = 0; i < m.rows (); i++)
239  accj.accum (m(i, j));
240 
241  res.xelem (j) = accj;
242  }
243 }
244 
245 template <class T, class R, class ACC>
246 void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
247 {
248  res = MArray<R> (dim_vector (m.rows (), 1));
249  std::vector<ACC> acci (m.rows (), acc);
250  for (octave_idx_type j = 0; j < m.columns (); j++)
251  {
252  for (octave_idx_type i = 0; i < m.rows (); i++)
253  acci[i].accum (m(i, j));
254  }
255 
256  for (octave_idx_type i = 0; i < m.rows (); i++)
257  res.xelem (i) = acci[i];
258 }
259 
260 // sparse versions
261 template <class T, class R, class ACC>
262 void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
263 {
264  res = MArray<R> (dim_vector (1, m.columns ()));
265  for (octave_idx_type j = 0; j < m.columns (); j++)
266  {
267  ACC accj = acc;
268  for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
269  accj.accum (m.data (k));
270 
271  res.xelem (j) = accj;
272  }
273 }
274 
275 template <class T, class R, class ACC>
276 void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
277 {
278  res = MArray<R> (dim_vector (m.rows (), 1));
279  std::vector<ACC> acci (m.rows (), acc);
280  for (octave_idx_type j = 0; j < m.columns (); j++)
281  {
282  for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
283  acci[m.ridx (k)].accum (m.data (k));
284  }
285 
286  for (octave_idx_type i = 0; i < m.rows (); i++)
287  res.xelem (i) = acci[i];
288 }
289 
290 // now the dispatchers
291 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \
292 template <class T, class R> \
293 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \
294 { \
295  RES_TYPE res; \
296  if (p == 2) \
297  FUNC_NAME (v, res, norm_accumulator_2<R> ()); \
298  else if (p == 1) \
299  FUNC_NAME (v, res, norm_accumulator_1<R> ()); \
300  else if (lo_ieee_isinf (p)) \
301  { \
302  if (p > 0) \
303  FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \
304  else \
305  FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \
306  } \
307  else if (p == 0) \
308  FUNC_NAME (v, res, norm_accumulator_0<R> ()); \
309  else if (p > 0) \
310  FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \
311  else \
312  FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \
313  return res; \
314 }
315 
319 DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>)
320 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>)
321 
322 // The approximate subproblem in Higham's method. Find lambda and mu such that
323 // norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized.
324 // Real version. As in Higham's paper.
325 template <class ColVectorT, class R>
326 static void
327 higham_subp (const ColVectorT& y, const ColVectorT& col,
328  octave_idx_type nsamp, R p, R& lambda, R& mu)
329 {
330  R nrm = 0;
331  for (octave_idx_type i = 0; i < nsamp; i++)
332  {
333  octave_quit ();
334  R fi = i * static_cast<R> (M_PI) / nsamp;
335  R lambda1 = cos (fi);
336  R mu1 = sin (fi);
337  R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
338  std::pow (std::abs (mu1), p), 1/p);
339  lambda1 /= lmnr; mu1 /= lmnr;
340  R nrm1 = vector_norm (lambda1 * y + mu1 * col, p);
341  if (nrm1 > nrm)
342  {
343  lambda = lambda1;
344  mu = mu1;
345  nrm = nrm1;
346  }
347  }
348 }
349 
350 // Complex version. Higham's paper does not deal with complex case, so we use a
351 // simple extension. First, guess the magnitudes as in real version, then try
352 // to rotate lambda to improve further.
353 template <class ColVectorT, class R>
354 static void
355 higham_subp (const ColVectorT& y, const ColVectorT& col,
356  octave_idx_type nsamp, R p,
357  std::complex<R>& lambda, std::complex<R>& mu)
358 {
359  typedef std::complex<R> CR;
360  R nrm = 0;
361  lambda = 1.0;
362  CR lamcu = lambda / std::abs (lambda);
363  // Probe magnitudes
364  for (octave_idx_type i = 0; i < nsamp; i++)
365  {
366  octave_quit ();
367  R fi = i * static_cast<R> (M_PI) / nsamp;
368  R lambda1 = cos (fi);
369  R mu1 = sin (fi);
370  R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
371  std::pow (std::abs (mu1), p), 1/p);
372  lambda1 /= lmnr; mu1 /= lmnr;
373  R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p);
374  if (nrm1 > nrm)
375  {
376  lambda = lambda1 * lamcu;
377  mu = mu1;
378  nrm = nrm1;
379  }
380  }
381  R lama = std::abs (lambda);
382  // Probe orientation
383  for (octave_idx_type i = 0; i < nsamp; i++)
384  {
385  octave_quit ();
386  R fi = i * static_cast<R> (M_PI) / nsamp;
387  lamcu = CR (cos (fi), sin (fi));
388  R nrm1 = vector_norm (lama * lamcu * y + mu * col, p);
389  if (nrm1 > nrm)
390  {
391  lambda = lama * lamcu;
392  nrm = nrm1;
393  }
394  }
395 }
396 
397 // the p-dual element (should work for both real and complex)
398 template <class T, class R>
399 inline T elem_dual_p (T x, R p)
400 {
401  return signum (x) * std::pow (std::abs (x), p-1);
402 }
403 
404 // the VectorT is used for vectors, but actually it has to be
405 // a Matrix type to allow all the operations. For instance SparseMatrix
406 // does not support multiplication with column/row vectors.
407 // the dual vector
408 template <class VectorT, class R>
409 VectorT dual_p (const VectorT& x, R p, R q)
410 {
411  VectorT res (x.dims ());
412  for (octave_idx_type i = 0; i < x.numel (); i++)
413  res.xelem (i) = elem_dual_p (x(i), p);
414  return res / vector_norm (res, q);
415 }
416 
417 // Higham's hybrid method
418 template <class MatrixT, class VectorT, class R>
419 R higham (const MatrixT& m, R p, R tol, int maxiter,
420  VectorT& x)
421 {
422  x.resize (m.columns (), 1);
423  // the OSE part
424  VectorT y(m.rows (), 1, 0), z(m.rows (), 1);
425  typedef typename VectorT::element_type RR;
426  RR lambda = 0;
427  RR mu = 1;
428  for (octave_idx_type k = 0; k < m.columns (); k++)
429  {
430  octave_quit ();
431  VectorT col (m.column (k));
432  if (k > 0)
433  higham_subp (y, col, 4*k, p, lambda, mu);
434  for (octave_idx_type i = 0; i < k; i++)
435  x(i) *= lambda;
436  x(k) = mu;
437  y = lambda * y + mu * col;
438  }
439 
440  // the PM part
441  x = x / vector_norm (x, p);
442  R q = p/(p-1);
443 
444  R gamma = 0, gamma1;
445  int iter = 0;
446  while (iter < maxiter)
447  {
448  octave_quit ();
449  y = m*x;
450  gamma1 = gamma;
451  gamma = vector_norm (y, p);
452  z = dual_p (y, p, q);
453  z = z.hermitian ();
454  z = z * m;
455 
456  if (iter > 0 && (vector_norm (z, q) <= gamma
457  || (gamma - gamma1) <= tol*gamma))
458  break;
459 
460  z = z.hermitian ();
461  x = dual_p (z, q, p);
462  iter ++;
463  }
464 
465  return gamma;
466 }
467 
468 // derive column vector and SVD types
469 
470 static const char *p_less1_gripe = "xnorm: p must be at least 1";
471 
472 // Static constant to control the maximum number of iterations. 100 seems to
473 // be a good value. Eventually, we can provide a means to change this
474 // constant from Octave.
475 static int max_norm_iter = 100;
476 
477 // version with SVD for dense matrices
478 template <class MatrixT, class VectorT, class SVDT, class R>
479 R matrix_norm (const MatrixT& m, R p, VectorT, SVDT)
480 {
481  R res = 0;
482  if (p == 2)
483  {
484  SVDT svd (m, SVD::sigma_only);
485  res = svd.singular_values () (0,0);
486  }
487  else if (p == 1)
488  res = xcolnorms (m, 1).max ();
489  else if (lo_ieee_isinf (p))
490  res = xrownorms (m, 1).max ();
491  else if (p > 1)
492  {
493  VectorT x;
494  const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
495  res = higham (m, p, sqrteps, max_norm_iter, x);
496  }
497  else
499 
500  return res;
501 }
502 
503 // SVD-free version for sparse matrices
504 template <class MatrixT, class VectorT, class R>
505 R matrix_norm (const MatrixT& m, R p, VectorT)
506 {
507  R res = 0;
508  if (p == 1)
509  res = xcolnorms (m, 1).max ();
510  else if (lo_ieee_isinf (p))
511  res = xrownorms (m, 1).max ();
512  else if (p > 1)
513  {
514  VectorT x;
515  const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
516  res = higham (m, p, sqrteps, max_norm_iter, x);
517  }
518  else
520 
521  return res;
522 }
523 
524 // and finally, here's what we've promised in the header file
525 
526 #define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \
527  OCTAVE_API RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \
528  { return vector_norm (x, p); } \
529  OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \
530  { return vector_norm (x, p); } \
531  OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \
532  { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \
533  OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \
534  { return vector_norm (x, static_cast<RTYPE> (2)); }
535 
538 DEFINE_XNORM_FUNCS(Float, float)
540 
541 // this is needed to avoid copying the sparse matrix for xfrobnorm
542 template <class T, class R>
543 inline void array_norm_2 (const T* v, octave_idx_type n, R& res)
544 {
546  for (octave_idx_type i = 0; i < n; i++)
547  acc.accum (v[i]);
548 
549  res = acc;
550 }
551 
552 #define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \
553  OCTAVE_API RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \
554  { return matrix_norm (x, p, PREFIX##Matrix ()); } \
555  OCTAVE_API RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \
556  { \
557  RTYPE res; \
558  array_norm_2 (x.data (), x.nnz (), res); \
559  return res; \
560  }
561 
564 
565 #define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \
566  extern OCTAVE_API RPREFIX##RowVector xcolnorms (const PREFIX##Matrix& m, RTYPE p) \
567  { return column_norms (m, p); } \
568  extern OCTAVE_API RPREFIX##ColumnVector xrownorms (const PREFIX##Matrix& m, RTYPE p) \
569  { return row_norms (m, p); } \
570 
572 DEFINE_COLROW_NORM_FUNCS(Complex, , double)
573 DEFINE_COLROW_NORM_FUNCS(Float, Float, float)
575 
577 DEFINE_COLROW_NORM_FUNCS(SparseComplex, , double)
578 
OCTAVE_API ColumnVector xrownorms(const Matrix &m, double p)
Definition: oct-norm.cc:571
#define DEFINE_XNORM_FUNCS(PREFIX, RTYPE)
Definition: oct-norm.cc:526
T * data(void)
Definition: Sparse.h:509
octave_idx_type rows(void) const
Definition: Sparse.h:263
#define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE)
Definition: oct-norm.cc:291
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:275
VectorT dual_p(const VectorT &x, R p, R q)
Definition: oct-norm.cc:409
void accum(U val)
Definition: oct-norm.cc:196
void accum(R val)
Definition: oct-norm.cc:134
octave_idx_type * cidx(void)
Definition: Sparse.h:531
octave_idx_type columns(void) const
Definition: Sparse.h:265
#define lo_ieee_isinf(x)
Definition: lo-ieee.h:105
Definition: MArray.h:36
R matrix_norm(const MatrixT &m, R p, VectorT, SVDT)
Definition: oct-norm.cc:479
void accum(U val)
Definition: oct-norm.cc:107
octave_idx_type rows(void) const
Definition: Array.h:313
T elem_dual_p(T x, R p)
Definition: oct-norm.cc:399
double signum(double x)
Definition: lo-mappers.cc:80
liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
static void higham_subp(const ColVectorT &y, const ColVectorT &col, octave_idx_type nsamp, R p, R &lambda, R &mu)
Definition: oct-norm.cc:327
void accum(U val)
Definition: oct-norm.cc:166
static int max_norm_iter
Definition: oct-norm.cc:475
#define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE)
Definition: oct-norm.cc:565
#define octave_Inf
Definition: lo-ieee.h:31
double max(void) const
Definition: dRowVector.cc:246
OCTAVE_API RowVector xcolnorms(const Matrix &m, double p)
Definition: oct-norm.cc:571
static R pow2(R x)
Definition: oct-norm.cc:130
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
norm_accumulator_p(R pp)
Definition: oct-norm.cc:76
unsigned int num
Definition: oct-norm.cc:207
static const char * p_less1_gripe
Definition: oct-norm.cc:470
void accum(std::complex< R > val)
Definition: oct-norm.cc:149
void accum(U val)
Definition: oct-norm.cc:211
T & xelem(octave_idx_type n)
Definition: Array.h:353
Handles the reference counting for all the derived classes.
Definition: Array.h:45
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
octave_idx_type * ridx(void)
Definition: Sparse.h:518
double max(void) const
Definition: dColVector.cc:284
#define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE)
Definition: oct-norm.cc:552
norm_accumulator_mp(R pp)
Definition: oct-norm.cc:104
void array_norm_2(const T *v, octave_idx_type n, R &res)
Definition: oct-norm.cc:543
static double fi[256]
Definition: randmtzig.c:443
function gamma(X)
Definition: gamma.f:2
void accum(U val)
Definition: oct-norm.cc:79
void row_norms(const MArray< T > &m, MArray< R > &res, ACC acc)
Definition: oct-norm.cc:246
void vector_norm(const Array< T > &v, R &res, ACC acc)
Definition: oct-norm.cc:222
void accum(U val)
Definition: oct-norm.cc:181
void column_norms(const MArray< T > &m, MArray< R > &res, ACC acc)
Definition: oct-norm.cc:232
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
R higham(const MatrixT &m, R p, R tol, int maxiter, VectorT &x)
Definition: oct-norm.cc:419
std::complex< double > Complex
Definition: oct-cmplx.h:29
T abs(T x)
Definition: pr-output.cc:3062
octave_idx_type columns(void) const
Definition: Array.h:322
F77_RET_T const double * x
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:210