GNU Octave  3.8.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-2013 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 <[email protected]>
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 "Array.h"
40 #include "Array-util.h"
41 #include "CMatrix.h"
42 #include "dMatrix.h"
43 #include "fCMatrix.h"
44 #include "fMatrix.h"
45 #include "CColVector.h"
46 #include "dColVector.h"
47 #include "CRowVector.h"
48 #include "dRowVector.h"
49 #include "fCColVector.h"
50 #include "fColVector.h"
51 #include "fCRowVector.h"
52 #include "fRowVector.h"
53 #include "CSparse.h"
54 #include "dSparse.h"
55 #include "dbleSVD.h"
56 #include "CmplxSVD.h"
57 #include "floatSVD.h"
58 #include "fCmplxSVD.h"
59 
60 // Theory: norm accumulator is an object that has an accum method able
61 // to handle both real and complex element, and a cast operator
62 // returning the intermediate norm. Reference: Higham, N. "Estimating
63 // the Matrix p-Norm." Numer. Math. 62, 539-555, 1992.
64 
65 // norm accumulator for the p-norm
66 template <class R>
68 {
69  R p,scl,sum;
70 public:
71  norm_accumulator_p () {} // we need this one for Array
72  norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {}
73 
74  template<class U>
75  void accum (U val)
76  {
77  octave_quit ();
78  R t = std::abs (val);
79  if (scl == t) // we need this to handle Infs properly
80  sum += 1;
81  else if (scl < t)
82  {
83  sum *= std::pow (scl/t, p);
84  sum += 1;
85  scl = t;
86  }
87  else if (t != 0)
88  sum += std::pow (t/scl, p);
89  }
90  operator R () { return scl * std::pow (sum, 1/p); }
91 };
92 
93 // norm accumulator for the minus p-pseudonorm
94 template <class R>
96 {
97  R p,scl,sum;
98 public:
99  norm_accumulator_mp () {} // we need this one for Array
100  norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {}
101 
102  template<class U>
103  void accum (U val)
104  {
105  octave_quit ();
106  R t = 1 / std::abs (val);
107  if (scl == t)
108  sum += 1;
109  else if (scl < t)
110  {
111  sum *= std::pow (scl/t, p);
112  sum += 1;
113  scl = t;
114  }
115  else if (t != 0)
116  sum += std::pow (t/scl, p);
117  }
118  operator R () { return scl * std::pow (sum, -1/p); }
119 };
120 
121 // norm accumulator for the 2-norm (euclidean)
122 template <class R>
124 {
125  R scl,sum;
126  static R pow2 (R x) { return x*x; }
127 public:
128  norm_accumulator_2 () : scl(0), sum(1) {}
129 
130  void accum (R val)
131  {
132  R t = std::abs (val);
133  if (scl == t)
134  sum += 1;
135  else if (scl < t)
136  {
137  sum *= pow2 (scl/t);
138  sum += 1;
139  scl = t;
140  }
141  else if (t != 0)
142  sum += pow2 (t/scl);
143  }
144 
145  void accum (std::complex<R> val)
146  {
147  accum (val.real ());
148  accum (val.imag ());
149  }
150 
151  operator R () { return scl * std::sqrt (sum); }
152 };
153 
154 // norm accumulator for the 1-norm (city metric)
155 template <class R>
157 {
158  R sum;
159 public:
161  template<class U>
162  void accum (U val)
163  {
164  sum += std::abs (val);
165  }
166  operator R () { return sum; }
167 };
168 
169 // norm accumulator for the inf-norm (max metric)
170 template <class R>
172 {
173  R max;
174 public:
176  template<class U>
177  void accum (U val)
178  {
179  max = std::max (max, std::abs (val));
180  }
181  operator R () { return max; }
182 };
183 
184 // norm accumulator for the -inf pseudonorm (min abs value)
185 template <class R>
187 {
188  R min;
189 public:
191  template<class U>
192  void accum (U val)
193  {
194  min = std::min (min, std::abs (val));
195  }
196  operator R () { return min; }
197 };
198 
199 // norm accumulator for the 0-pseudonorm (hamming distance)
200 template <class R>
202 {
203  unsigned int num;
204 public:
206  template<class U>
207  void accum (U val)
208  {
209  if (val != static_cast<U> (0)) ++num;
210  }
211  operator R () { return num; }
212 };
213 
214 
215 // OK, we're armed :) Now let's go for the fun
216 
217 template <class T, class R, class ACC>
218 inline void vector_norm (const Array<T>& v, R& res, ACC acc)
219 {
220  for (octave_idx_type i = 0; i < v.numel (); i++)
221  acc.accum (v(i));
222 
223  res = acc;
224 }
225 
226 // dense versions
227 template <class T, class R, class ACC>
228 void column_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
229 {
230  res = MArray<R> (dim_vector (1, m.columns ()));
231  for (octave_idx_type j = 0; j < m.columns (); j++)
232  {
233  ACC accj = acc;
234  for (octave_idx_type i = 0; i < m.rows (); i++)
235  accj.accum (m(i, j));
236 
237  res.xelem (j) = accj;
238  }
239 }
240 
241 template <class T, class R, class ACC>
242 void row_norms (const MArray<T>& m, MArray<R>& res, ACC acc)
243 {
244  res = MArray<R> (dim_vector (m.rows (), 1));
245  std::vector<ACC> acci (m.rows (), acc);
246  for (octave_idx_type j = 0; j < m.columns (); j++)
247  {
248  for (octave_idx_type i = 0; i < m.rows (); i++)
249  acci[i].accum (m(i, j));
250  }
251 
252  for (octave_idx_type i = 0; i < m.rows (); i++)
253  res.xelem (i) = acci[i];
254 }
255 
256 // sparse versions
257 template <class T, class R, class ACC>
258 void column_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
259 {
260  res = MArray<R> (dim_vector (1, m.columns ()));
261  for (octave_idx_type j = 0; j < m.columns (); j++)
262  {
263  ACC accj = acc;
264  for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
265  accj.accum (m.data (k));
266 
267  res.xelem (j) = accj;
268  }
269 }
270 
271 template <class T, class R, class ACC>
272 void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc)
273 {
274  res = MArray<R> (dim_vector (m.rows (), 1));
275  std::vector<ACC> acci (m.rows (), acc);
276  for (octave_idx_type j = 0; j < m.columns (); j++)
277  {
278  for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++)
279  acci[m.ridx (k)].accum (m.data (k));
280  }
281 
282  for (octave_idx_type i = 0; i < m.rows (); i++)
283  res.xelem (i) = acci[i];
284 }
285 
286 // now the dispatchers
287 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \
288 template <class T, class R> \
289 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \
290 { \
291  RES_TYPE res; \
292  if (p == 2) \
293  FUNC_NAME (v, res, norm_accumulator_2<R> ()); \
294  else if (p == 1) \
295  FUNC_NAME (v, res, norm_accumulator_1<R> ()); \
296  else if (lo_ieee_isinf (p)) \
297  { \
298  if (p > 0) \
299  FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \
300  else \
301  FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \
302  } \
303  else if (p == 0) \
304  FUNC_NAME (v, res, norm_accumulator_0<R> ()); \
305  else if (p > 0) \
306  FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \
307  else \
308  FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \
309  return res; \
310 }
311 
315 DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>)
316 DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>)
317 
318 // The approximate subproblem in Higham's method. Find lambda and mu such that
319 // norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized.
320 // Real version. As in Higham's paper.
321 template <class ColVectorT, class R>
322 static void
323 higham_subp (const ColVectorT& y, const ColVectorT& col,
324  octave_idx_type nsamp, R p, R& lambda, R& mu)
325 {
326  R nrm = 0;
327  for (octave_idx_type i = 0; i < nsamp; i++)
328  {
329  octave_quit ();
330  R fi = i*M_PI/nsamp, lambda1 = cos (fi), mu1 = sin (fi);
331  R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
332  std::pow (std::abs (mu1), p), 1/p);
333  lambda1 /= lmnr; mu1 /= lmnr;
334  R nrm1 = vector_norm (lambda1 * y + mu1 * col, p);
335  if (nrm1 > nrm)
336  {
337  lambda = lambda1;
338  mu = mu1;
339  nrm = nrm1;
340  }
341  }
342 }
343 
344 // Complex version. Higham's paper does not deal with complex case, so we use a
345 // simple extension. First, guess the magnitudes as in real version, then try
346 // to rotate lambda to improve further.
347 template <class ColVectorT, class R>
348 static void
349 higham_subp (const ColVectorT& y, const ColVectorT& col,
350  octave_idx_type nsamp, R p,
351  std::complex<R>& lambda, std::complex<R>& mu)
352 {
353  typedef std::complex<R> CR;
354  R nrm = 0;
355  lambda = 1.0;
356  CR lamcu = lambda / std::abs (lambda);
357  // Probe magnitudes
358  for (octave_idx_type i = 0; i < nsamp; i++)
359  {
360  octave_quit ();
361  R fi = i*M_PI/nsamp, lambda1 = cos (fi), mu1 = sin (fi);
362  R lmnr = std::pow (std::pow (std::abs (lambda1), p) +
363  std::pow (std::abs (mu1), p), 1/p);
364  lambda1 /= lmnr; mu1 /= lmnr;
365  R nrm1 = vector_norm (lambda1 * lamcu * y + mu1 * col, p);
366  if (nrm1 > nrm)
367  {
368  lambda = lambda1 * lamcu;
369  mu = mu1;
370  nrm = nrm1;
371  }
372  }
373  R lama = std::abs (lambda);
374  // Probe orientation
375  for (octave_idx_type i = 0; i < nsamp; i++)
376  {
377  octave_quit ();
378  R fi = i*M_PI/nsamp;
379  lamcu = CR (cos (fi), sin (fi));
380  R nrm1 = vector_norm (lama * lamcu * y + mu * col, p);
381  if (nrm1 > nrm)
382  {
383  lambda = lama * lamcu;
384  nrm = nrm1;
385  }
386  }
387 }
388 
389 // the p-dual element (should work for both real and complex)
390 template <class T, class R>
391 inline T elem_dual_p (T x, R p)
392 {
393  return signum (x) * std::pow (std::abs (x), p-1);
394 }
395 
396 // the VectorT is used for vectors, but actually it has to be
397 // a Matrix type to allow all the operations. For instance SparseMatrix
398 // does not support multiplication with column/row vectors.
399 // the dual vector
400 template <class VectorT, class R>
401 VectorT dual_p (const VectorT& x, R p, R q)
402 {
403  VectorT res (x.dims ());
404  for (octave_idx_type i = 0; i < x.numel (); i++)
405  res.xelem (i) = elem_dual_p (x(i), p);
406  return res / vector_norm (res, q);
407 }
408 
409 // Higham's hybrid method
410 template <class MatrixT, class VectorT, class R>
411 R higham (const MatrixT& m, R p, R tol, int maxiter,
412  VectorT& x)
413 {
414  x.resize (m.columns (), 1);
415  // the OSE part
416  VectorT y(m.rows (), 1, 0), z(m.rows (), 1);
417  typedef typename VectorT::element_type RR;
418  RR lambda = 0, mu = 1;
419  for (octave_idx_type k = 0; k < m.columns (); k++)
420  {
421  octave_quit ();
422  VectorT col (m.column (k));
423  if (k > 0)
424  higham_subp (y, col, 4*k, p, lambda, mu);
425  for (octave_idx_type i = 0; i < k; i++)
426  x(i) *= lambda;
427  x(k) = mu;
428  y = lambda * y + mu * col;
429  }
430 
431  // the PM part
432  x = x / vector_norm (x, p);
433  R q = p/(p-1);
434 
435  R gamma = 0, gamma1;
436  int iter = 0;
437  while (iter < maxiter)
438  {
439  octave_quit ();
440  y = m*x;
441  gamma1 = gamma;
442  gamma = vector_norm (y, p);
443  z = dual_p (y, p, q);
444  z = z.hermitian ();
445  z = z * m;
446 
447  if (iter > 0 && (vector_norm (z, q) <= gamma
448  || (gamma - gamma1) <= tol*gamma))
449  break;
450 
451  z = z.hermitian ();
452  x = dual_p (z, q, p);
453  iter ++;
454  }
455 
456  return gamma;
457 }
458 
459 // derive column vector and SVD types
460 
461 static const char *p_less1_gripe = "xnorm: p must be at least 1";
462 
463 // Static constant to control the maximum number of iterations. 100 seems to
464 // be a good value. Eventually, we can provide a means to change this
465 // constant from Octave.
466 static int max_norm_iter = 100;
467 
468 // version with SVD for dense matrices
469 template <class MatrixT, class VectorT, class SVDT, class R>
470 R matrix_norm (const MatrixT& m, R p, VectorT, SVDT)
471 {
472  R res = 0;
473  if (p == 2)
474  {
475  SVDT svd (m, SVD::sigma_only);
476  res = svd.singular_values () (0,0);
477  }
478  else if (p == 1)
479  res = xcolnorms (m, 1).max ();
480  else if (lo_ieee_isinf (p))
481  res = xrownorms (m, 1).max ();
482  else if (p > 1)
483  {
484  VectorT x;
485  const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
486  res = higham (m, p, sqrteps, max_norm_iter, x);
487  }
488  else
490 
491  return res;
492 }
493 
494 // SVD-free version for sparse matrices
495 template <class MatrixT, class VectorT, class R>
496 R matrix_norm (const MatrixT& m, R p, VectorT)
497 {
498  R res = 0;
499  if (p == 1)
500  res = xcolnorms (m, 1).max ();
501  else if (lo_ieee_isinf (p))
502  res = xrownorms (m, 1).max ();
503  else if (p > 1)
504  {
505  VectorT x;
506  const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
507  res = higham (m, p, sqrteps, max_norm_iter, x);
508  }
509  else
511 
512  return res;
513 }
514 
515 // and finally, here's what we've promised in the header file
516 
517 #define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \
518  OCTAVE_API RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \
519  { return vector_norm (x, p); } \
520  OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \
521  { return vector_norm (x, p); } \
522  OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \
523  { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \
524  OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \
525  { return vector_norm (x, static_cast<RTYPE> (2)); }
526 
529 DEFINE_XNORM_FUNCS(Float, float)
531 
532 // this is needed to avoid copying the sparse matrix for xfrobnorm
533 template <class T, class R>
534 inline void array_norm_2 (const T* v, octave_idx_type n, R& res)
535 {
537  for (octave_idx_type i = 0; i < n; i++)
538  acc.accum (v[i]);
539 
540  res = acc;
541 }
542 
543 #define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \
544  OCTAVE_API RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \
545  { return matrix_norm (x, p, PREFIX##Matrix ()); } \
546  OCTAVE_API RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \
547  { \
548  RTYPE res; \
549  array_norm_2 (x.data (), x.nnz (), res); \
550  return res; \
551  }
552 
555 
556 #define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \
557  extern OCTAVE_API RPREFIX##RowVector xcolnorms (const PREFIX##Matrix& m, RTYPE p) \
558  { return column_norms (m, p); } \
559  extern OCTAVE_API RPREFIX##ColumnVector xrownorms (const PREFIX##Matrix& m, RTYPE p) \
560  { return row_norms (m, p); } \
561 
563 DEFINE_COLROW_NORM_FUNCS(Complex, , double)
564 DEFINE_COLROW_NORM_FUNCS(Float, Float, float)
566 
568 DEFINE_COLROW_NORM_FUNCS(SparseComplex, , double)
569