From octave-sources-request at bevo dot che dot wisc dot edu  Tue Oct  3 07:06:55 2000
Subject: Re: function : outer product
From: Paul Kienzle <pkienzle at kienzle dot powernet dot co dot uk>
To: Antoine Mathys <mathyant at student dot ethz dot ch>
Cc: octave-sources at bevo dot che dot wisc dot edu
Date: Tue, 3 Oct 2000 12:32:02 +0100

feval(f, x(i), y(j)) is about 2.5 times faster than eval(sprintf(...)),
and it doesn't cause a conflict with local variables (I happened to
choose 'f' as my function, which in the eval context evaluates to the
string 'f' rather than the function f().

The direct outer product can be more much more quickly computed, though
with three times the memory, by building [x,x,x,...] and [y;y;y;...]
and doing element-wise multiplication like so:

    x = x (:);
    x = x (:, ones (1, n))
    y = y (:).';
    y = y (ones (1, m), :)
    retval = x .* y;

If we know that the function could accept matrix arguments (perhaps with
an extra parameter?), we can use the same technique.

Note that the memory problem could be removed if octave supported M .*=
v, allowing us to write:

    retval = ones (m, n);
    retval .*= x(:);
    retval .*= y(:).';

Paul Kienzle
pkienzle at kienzle dot powernet dot co dot uk

## Copyright (C) 2000 Antoine Mathys
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## -*- texinfo -*-
##  at deftypefn {Function File} {@var{a} =} outer (@var{x}, @var{y}, [@var{f}])
## Given two vectors of length  at var{m} and @var{n}, return a
##  at var{m}-by-@var{n} matrix with entries @var{f (x(i), y(j))}.
##
## If  at var{f} is omitted, multiplication is the default.
##
##  at deftypefnx {Function File} {@var{a} =} outer (@var{x}, @var{y}, @var{f}, 'matrix')
## Pass matrices to function  at var{f} instead of individual elements.
##  at end deftypefn

## Author: Antoine Mathys <mathyant at student dot ethz dot ch>
## Created: 2 October 2000
## Adapted-By: Antoine Mathys

function retval = outer (x, y, f, matrix)

  if (nargin < 2 || nargin > 4)
    usage ("outer (x, y [, f [, 'matrix']])");
  endif
  
  if (!is_vector(x) || !is_vector(y))
    error ("outer: both x and y must be vectors");
  endif

  m = length (x);
  n = length (y);

  if nargin == 2
    x = x (:);
    x = x (:, ones (1, n));
    y = y (:).';
    y = y (ones (1, m), :);
    retval = x .* y;
  elseif nargin == 4
    x = x (:);
    x = x (:, ones (1, n));
    y = y (:).';
    y = y (ones (1, m), :);
    retval = feval (f, x, y);
  else
    retval = zeros (m, n);
    for i=1:m
      for j=1:n
	retval (i, j) = feval (f, x(i), y(j));
      endfor
    endfor
  endif

endfunction


On Mon, Oct 02, 2000 at 06:47:15PM +0200, Antoine Mathys wrote:
> Hello,
> 
> Here is an outer function, as described in octave's project file.
> outer(x,y,f) returns a matrix with f(x(i),y(j)) as elements.
> 
> 
> 
> --
> --------------------------------------------------
> Antoine Mathys / e-mail : mathyant at student dot ethz dot ch
> URL : http://n.ethz.ch/student/mathyant/
> --------------------------------------------------
> 
> 

> ## Copyright (C) 2000 Antoine Mathys
> ##
> ## This file is part of Octave.
> ##
> ## Octave is free software; you can redistribute it and/or modify it
> ## under the terms of the GNU General Public License as published by
> ## the Free Software Foundation; either version 2, or (at your option)
> ## any later version.
> ##
> ## Octave is distributed in the hope that it will be useful, but
> ## WITHOUT ANY WARRANTY; without even the implied warranty of
> ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> ## General Public License for more details.
> ##
> ## You should have received a copy of the GNU General Public License
> ## along with Octave; see the file COPYING.  If not, write to the Free
> ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
> ## 02111-1307, USA.
> 
> ## -*- texinfo -*-
> ##  at deftypefn {Function File} {@var{a} =} outer (@var{x}, @var{y}, [@var{f}])
> ## Given two vectors of length  at var{m} and @var{n}, return a
> ##  at var{m}-by-@var{n} matrix with entries @var{f (x(i), y(j))}.
> ##
> ## If  at var{f} is omitted, multiplication is the default.
> ##  at end deftypefn
> 
> ## Author: Antoine Mathys <mathyant at student dot ethz dot ch>
> ## Created: 2 October 2000
> ## Adapted-By: Antoine Mathys
> 
> function retval = outer (x, y, f)
> 
>   if (nargin < 2 || nargin > 3)
>     usage ("outer (x, y, [f])");
>   endif
>   
>   if (!is_vector(x) || !is_vector(y))
>     error ("outer: both x and y must be vectors");
>   endif
> 
>   m = length (x);
>   n = length (y);
> 
>   retval = zeros (m, n);
>   
>   for i=1:m
>     for j=1:n
>       if (nargin == 2)
> 	retval(i,j) = x(i) * y(j);
>       else
> 	retval(i,j) = eval(sprintf("%s (x(i), y(j));",f));
>       endif
>     endfor
>   endfor
> 
> endfunction



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------