# [SciPy-user] Numerical gradient approximation on matrix

Travis Oliphant oliphant at ee.byu.edu
Mon Aug 1 15:02:59 CDT 2005

```Dimitri D'Or wrote:

>> Message: 7
>
>> Date: Fri, 29 Jul 2005 14:42:12 -0400
>
>> From: Alan G Isaac <aisaac at american.edu>
>
>> Subject: Re[2]: [SciPy-user] Numerical gradient approximation on
>
>>     matrix
>
>> To: SciPy Users List <scipy-user at scipy.net>
>
>> Message-ID: <Mahogany-0.66.0-1428-20050729-144212.01 at american.edu>
>
>> Content-Type: TEXT/PLAIN; CHARSET=UTF-8
>
>>
>
>> On Fri, 29 Jul 2005, T) guillem at torroja.dmt.upm.es apparently wrote:
>
>
>> > and keeps the array's shape:
>
>>
>
>> How's that supposed to work.
>
>> I'm not a Matlab user,
>
>> but the docs say it uses diff.
>
>>
>
>> Cheers,
>
>> Alan Isaac
>
>
>
> The Matlab gradient function doesn't use diff. It raises similar but
> not identical results to diff. A short example to illustrate this:
>
>
>
Attached is a MATLAB-compatible N-dimensional gradient function using
central differences on the interior and forward/backward differences on
the edges.

Notice the use of slice objects to generalize to N-dimensions.  This
technique can be used in many situations and is quite handy. It can make
the code a little harder to read.

-Travis

-------------- next part --------------

N = len(f.shape)  # number of dimensions
n = len(varargs)
if n==0:
dx = [1.0]*N
elif n==1:
dx = [varargs[0]]*N
elif n==N:
dx = list(varargs)
else:
raise SyntaxError, "Invalid number of arguments"

# use central differences on interior and first differences on endpoints

outvals = []

# create slice objects --- initially all are [:,:,...,:]
slices = [None]*3
for k in range(3):
slices[k] = []
for j in range(N):
slices[k].append(slice(None))

for axis in range(N):
# select out appropriate parts for this dimension
out = zeros(f.shape, f.typecode())
slices[0][axis] = slice(1,-1)
slices[1][axis] = slice(2,None)
slices[2][axis] = slice(None,-2)
# 1d equivalent -- out[1:-1] = (f[2:] - f[:-2])/2.0
out[slices[0]] = (f[slices[1]] - f[slices[2]])/2.0
slices[0][axis] = 0
slices[1][axis] = 1
slices[2][axis] = 0
# 1d equivalent -- out[0] = (f[1] - f[0])
out[slices[0]] = (f[slices[1]] - f[slices[2]])
slices[0][axis] = -1
slices[1][axis] = -1
slices[2][axis] = -2
# 1d equivalent -- out[-1] = (f[-1] - f[-2])
out[slices[0]] = (f[slices[1]] - f[slices[2]])

# divide by step function
outvals.append(out / dx[axis])

# reset the slice object in this dimension to ":"
slices[0][axis] = slice(None)
slices[1][axis] = slice(None)
slices[2][axis] = slice(None)

if N == 1:
return outvals[0]
else:
return outvals

```