octave_packages/m/sparse/sprandsym.m

   1 ## Copyright (C) 2004-2012 David Bateman and Andy Adler
   2 ## Copyright (C) 2012 Jordi Gutiérrez Hermoso
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {} sprandsym (@var{n}, @var{d})
  22 ## @deftypefnx {Function File} {} sprandsym (@var{s})
  23 ## Generate a symmetric random sparse matrix.  The size of the matrix will be
  24 ## @var{n} by @var{n}, with a density of values given by @var{d}.
  25 ## @var{d} should be between 0 and 1. Values will be normally
  26 ## distributed with mean of zero and variance 1.
  27 ##
  28 ## If called with a single matrix argument, a random sparse matrix is
  29 ## generated wherever the matrix @var{S} is non-zero in its lower
  30 ## triangular part.
  31 ## @seealso{sprand, sprandn}
  32 ## @end deftypefn
  33
  34 function S = sprandsym (n, d)
  35
  36   if (nargin != 1 && nargin != 2)
  37     print_usage ();
  38   endif
  39
  40   if (nargin == 1)
  41     [i, j] = find (tril (n));
  42     [nr, nc] = size (n);
  43     S = sparse (i, j, randn (size (i)), nr, nc);
  44     S = S + tril (S, -1)';
  45     return;
  46   endif
  47
  48   if (!(isscalar (n) && n == fix (n) && n > 0))
  49     error ("sprandsym: N must be an integer greater than 0");
  50   endif
  51
  52   if (d < 0 || d > 1)
  53     error ("sprandsym: density D must be between 0 and 1");
  54   endif
  55
  56   ## Actual number of nonzero entries
  57   k = round (n^2*d);
  58
  59   ## Diagonal nonzero entries, same parity as k
  60   r = pick_rand_diag (n, k);
  61
  62   ## Off diagonal nonzero entries
  63   m = (k - r)/2;
  64
  65   ondiag = randperm (n, r);
  66   offdiag = randperm (n*(n - 1)/2, m);
  67
  68   ## Row index
  69   i = lookup (cumsum (0:n), offdiag - 1) + 1;
  70
  71   ## Column index
  72   j = offdiag - (i - 1).*(i - 2)/2;
  73
  74   diagvals = randn (1, r);
  75   offdiagvals = randn (1, m);
  76
  77   S = sparse ([ondiag, i, j], [ondiag, j, i],
  78               [diagvals, offdiagvals, offdiagvals], n, n);
  79
  80 endfunction
  81
  82 function r = pick_rand_diag (n, k)
  83   ## Pick a random number R of entries for the diagonal of a sparse NxN
  84   ## symmetric square matrix with exactly K nonzero entries, ensuring
  85   ## that this R is chosen uniformly over all such matrices.
  86   ##
  87   ## Let D be the number of diagonal entries and M the number of
  88   ## off-diagonal entries. Then K = D + 2*M. Let A = N*(N-1)/2 be the
  89   ## number of available entries in the upper triangle of the matrix.
  90   ## Then, by a simple counting argument, there is a total of
  91   ##
  92   ##     T = nchoosek (N, D) * nchoosek (A, M)
  93   ##
  94   ## symmetric NxN matrices with a total of K nonzero entries and D on
  95   ## the diagonal. Letting D range from mod (K,2) through min (N,K), and
  96   ## dividing by this sum, we obtain the probability P for D to be each
  97   ## of those values.
  98   ##
  99   ## However, we cannot use this form for computation, as the binomial
 100   ## coefficients become unmanageably large. Instead, we use the
 101   ## successive quotients Q(i) = T(i+1)/T(i), which we easily compute to
 102   ## be
 103   ##
 104   ##               (N - D)*(N - D - 1)*M
 105   ##     Q =  -------------------------------
 106   ##            (D + 2)*(D + 1)*(A - M + 1)
 107   ##
 108   ## Then, after prepending 1, the cumprod of these quotients is
 109   ##
 110   ##      C = [ T(1)/T(1), T(2)/T(1), T(3)/T(1), ..., T(N)/T(1) ]
 111   ##
 112   ## Their sum is thus S = sum (T)/T(1), and then C(i)/S is the desired
 113   ## probability P(i) for i=1:N. The cumsum will finally give the
 114   ## distribution function for computing the random number of entries on
 115   ## the diagonal R.
 116   ##
 117   ## Thanks to Zsbán Ambrus <ambrus@math.bme.hu> for most of the ideas
 118   ## of the implementation here, especially how to do the computation
 119   ## numerically to avoid overflow.
 120
 121   ## Degenerate case
 122   if k == 1
 123     r = 1;
 124     return
 125   endif
 126
 127   ## Compute the stuff described above
 128   a = n*(n - 1)/2;
 129   d = [mod(k,2):2:min(n,k)-2];
 130   m = (k - d)/2;
 131   q = (n - d).*(n - d - 1).*m ./ (d + 2)./(d + 1)./(a - m + 1);
 132
 133   ## Slight modification from discussion above: pivot around the max in
 134   ## order to avoid overflow (underflow is fine, just means effectively
 135   ## zero probabilities).
 136   [~, midx] = max (cumsum (log (q))) ;
 137   midx++;
 138   lc = fliplr (cumprod (1./q(midx-1:-1:1)));
 139   rc = cumprod (q(midx:end));
 140
 141   ## Now c = t(i)/t(midx), so c > 1 == [].
 142   c = [lc, 1, rc];
 143   s = sum (c);
 144   p = c/s;
 145
 146   ## Add final d
 147   d(end+1) = d(end) + 2;
 148
 149   ## Pick a random r using this distribution
 150   r = d(sum (cumsum (p) < rand) + 1);
 151
 152 endfunction
 153
 154 %!test
 155 %! s = sprandsym (10, 0.1);
 156 %! assert (issparse (s));
 157 %! assert (issymmetric (s));
 158 %! assert (size (s), [10, 10]);
 159 %! assert (nnz (s) / numel (s), 0.1, .01);
 160
 161 %% Test 1-input calling form
 162 %!test
 163 %! s = sprandsym (sparse ([1 2 3], [3 2 3], [2 2 2]));
 164 %! [i, j] = find (s);
 165 %! assert (sort (i), [2 3]');
 166 %! assert (sort (j), [2 3]');
 167
 168 %% Test input validation
 169 %!error sprandsym ()
 170 %!error sprandsym (1, 2, 3)
 171 %!error sprandsym (ones(3), 0.5)
 172 %!error sprandsym (3.5, 0.5)
 173 %!error sprandsym (0, 0.5)
 174 %!error sprandsym (3, -1)
 175 %!error sprandsym (3, 2)
 176