The Winograd procedure for matrix multiplication provides a method whereby general matrix products may be computed more efficiently than the normal method. The algorithm and the time savings that can be effected are described. A FORTRAN program is provided which performs a general matrix multiply according to this algorithm. A variation of this procedure that may be used to calculate Gaussian probability density functions is also described. It is shown how a time savings can be effected in this calculation. The extension of this method to other similar calculations should yield similar savings

Lynn, M. S.

Schafer, R. A.

Vanrooy, D. L.

English

NASA Technical Reports Server

General Disclaimer 
One or more of the Following Statements may affect this Document 
 
 This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 
 
 This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 
 
 This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 
 
 This document is paginated as submitted by the original source. 
 
 Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 
 
 
 
 
 
 
 
Produced by the NASA Center for Aerospace Information (CASI) 
https://ntrs.nasa.gov/search.jsp?R=19760020815 2020-03-22T14:58:11+00:00Z
^	
1
NASA GR
(N)SA-CF-147C23)
	 TEE ISE CF THE WINCGFAD
	
N76-27903
NAIPIX PULTIFLICATICN ILGC Q ITHNJ IN rIGITIL
MULTISPICTEFL FFCCE^SING (Bice Univ.) 	 18 F
HC S3.5C	 C^CL (SSE	 Unclas
	G3/6
	 44565
.SAO)"..
	
'	 JUL 1976
^CEIlu
INPUT Bp-
INSTITUTE FOR CO'APUTER SERVICES AND APPLICATIONS
RICE, -. UNIVERSITY
275-025-027
The Use of the Winograd Matrix Multiplication
Algorithm in Dipital Multi-spectral Processing
by
1). L. Van Roov, R. A. Schafer, and
M. S. Lynn
1 C S A
Rice University
ABSTRACT:
The Winograd procedure for matrix multiplication. [S. Winograd,
Comm. on Pure and Applied Math. , 23, 19701 provid.-s a method
whereby  general matrix products may be computed more efficiently
than the normal method. In this report, we describe the algorithm
and the time savin rs that can be effected. , FORTRAN program is
provided which peXrms a general matrix muitiply according to this
algorithm.
Additionall y , we describe a variation of this procedure that may
be used to calCUlate Gaussian probability density functions. It is
shown how a time savingThes c n be effected in this calculation. ^e
extension of this inethocl to ether similar calculations should yield
similar savings.
Institute for Computer Services and Applications
Rice University
Houston, TX 77001
May, 1976
TLis work was supported under NASA contract i4AS 9-12776
1.	 lilt rmlut U( )n :
In this paper, we show the Winograd(  procedure for
computing n ► atrix products can be applied to various calculations
used in digital processing of remotely sensed data. The basic
procedure is described, a I'URTRAN program for general matrix
multiplication is provided, and an exanil-!e (computation of Gaussian
}p robability density functions) is worked out showing; the regions
where it	 -Imputatio; ally faster and the amount of time savings
involved. Of course, there are many other calculations where this
procedure can be effectively utilized.
Essentially, the Winograd procedure effects a time savings in
computing matrix products by trading off some of the multiplications
involved for additions (multiplies are usually slower operations on a
computer than adds). A relatively small amount of additional storage
is required, but a significant decrease in the member of multiplies
(tip to a factor of 2 for the case of multiplying two n x ii matrices,
if	 n	 is Gufficiertly large) can be gained. However, one must
sacrifice some numerical stability in employing this procedure(2).
In the next section, we describe the general matrix multipli-
cation procedure developed by Winograd and compare it with the
standard method.	 In Section 111 ; we uescribe a variation of this
procedure for the computation of Gaussian probability density
functions. An appendix containing a FORTRAN progru m for general
matrix multiplication is provided.
2II. The Winograd Algorithm :
Let A =(a ij ) he a p x q matrJ%, B=(b il ), a g x s
matrix, C= (c i j ) , a p x s matrix, x= l x i ), a q-vector,
and y = (y i) 	 a p-vector. The standard method for computing
y = A x is
(I
Al i	 I	 a i	 x	 i-1,2,	 ...pl	 j
j=1
and for computing C=AB is
q
cij	 -- a i ^l b qj	 i=1,2,...,p	 (2)
=1	 j=1,2,...,s
Thus to compute y in this manner, p q multiplies and p (q - 1 )
adds are necessary; and to compute C, p q s multiplies and
p s (q - 1 ) adds,
The Winograd procedure consists of rewriting eqs. (1) and (2)
so that 4ome quantities are precomputed. 'Mis procedure is based
on the identity
a ik b 	 + a i, k+ 1 bk+1
= (a ik + b k+1 )(a i, k+1 + 1'k ) - aik a i, k+1 - 1'k bk+1
Similar identities of a Nigher order may be used to construct other
algorithms, but, for out purposes here, these• are not of much
interest. Following Winograd's notation, we let I d J = largest
1
1integer
	 s d and	 f d l= sma llest integer	 > d,	 with	 E=	 f,	 q J,
the Winograd I)roceduty , then,	 for conilming y
3
\e
u=1
u=1
ai, 2u- I	 '	 a i, 2u
X 2u-I ' x 2u
i=1,2,...,p	 (3a)
(3b)
and
I	 ^-(a i, 2u- 1	 + X 2u )(a i, 2ti	 + x 2u- I )
u=1
y i =	 ^i - r	 if c1=2 e
( a` i, 2u- I	 + X 2st )(a i, 2u + X 2u- I) -
u=1
^ i - Ti - a i q x 	 if q=21 i i	 (3c)
This algor°thm requires p I q 1 4- 1(p+I) multiplications and
p ( 1 1, q I + ,e- I adds.	 if we have t y 's to compute using the
same A and t x's, the operation counts become
i p --- t - (p I"L( I
 I + t ) multi}] ies and p ( A - I ) + t - (I - I A
p(21+1+1^yl)^
recomputed. Table 1
etc. ) time savings to
various values of p
tit-no
	 'I	 to the adi
n^
for an IBM 370/155).
adds, since the p i 's need not be
shows the approximate (ignoring indexing,
be expected for computing 100 y's for
anti q for a ratio of the machine multiply
I time T t of 2.7 (the approximate value
We see from this table that for low values
a,
C,
w
w,U.s0$"o
7, Czv
p u ^
v UO
M!1 J m
a v
>
U
O U `
a,	 E
E r^
ct
U CO
o^
o •^ 3
03
cc
G'
W
RS
r
N
rr
a
w
t.
U
C:O
3
cc
bk
Cr
v
cz
4
N,
00 w r` to r` ID
^p In N w -c to ^! ^!' cn
O CT CO W o0 OG 00 aD
-- O O O O O G O
M N Cr UM 01 n v)
r- p Cr x o0 oc a0 Oo 00
-- O O O O O O O
00 t- t` to t • r-.
N 10 M CT r- IC to
00 00 00 00 00 00
^ -^ o 0 0 0 o c o
M to CT :n
O 0000 0000 0000 0000O a 00
-- O G O O O O O
W ^ ^ U J t^ !•'I t^
00OC Q T CT 00 00 w 00
^ O O O O O O O
oG 00 00 10 00 M 00
C^ M --+ O CT CT 00
^G ^' CT Cr CT C1^ 00 00 00
o 0 o c o 0 0
.-^ N 110 tV
N OC -Y co M
-r O CT C> CT CT CT 0
..; ^ O O O O !J O
^D ^p N CT Cl 00 00
N N ^ .-. + O O O C
Q' N er "D 00
Q
of p and q , therc is a net loss '-"e for large values, up to
,a
	 savings results.
The extension Of this concept to tell matrix multiplication is
(for C = A B)
P
	
^' 1	 1, 2u- I	 1, 2u
U = 1	 (4a)
I
	
j	 J	 b t u - I , j L'2u,j
u=1	 (4b)
t
2u I + b 2u, j )(a i, 2u + b2u
U=1
r; i - ri j 	if q =2p
2u- I + b 2u, I )(a i, 2u + b 2 i _ 1	 ^) -
u =1
^i	
- T) j +	 1 i71 1)	 if q =21, +- 1
for	 i= I , 2, ... , p	 and	 j=1 , 2, ... , s	 (4c)
11iis requires p s ['(I I + t (p + s) multiplies and ( i -1) (p +s) +
p s([ j g I+
	 2+	 1)	 adds. For	 p= q= s= n	 and n	 large,
this reduces	 to	 — A n 3 + n 2	 multiplies and
3
2 n 3 + n 2	 ides,
which can effect a savings in computation time over the standard
procedure which requires 	 n 3 multiplies and n 3 - n 2 adds;
using the ratio of 2.7: 1 for multiplies to adds, in fact a savings
always occurs for n > 5 .
6A FOR'T'RAN subroutine WNOMUL which performs matrix
multiplication according to this algorithm is listed in the appendix.
s
i[11. Calculation of Gaussian Probability Density functions :
The n-dimensional probability de .ity for a normnl population
with mean rl and covariance K is given by
n	 -
f(x) = (2-) W. I K	 exh I(x-u) T K - I (x-u) J
In this section, we shall concern ourselves only with the calculation
of the quadratic form (the argument of exp), given u and 1, and
D where L and 1) are the modified Cholesky decomposition(3)
of K, i. e.
	 K = L DL 3 ,
 
with L being unit lower triangular
and D, diagonal with positive diagonal entries. Then we can write
T
( x -11) T K - I (X-LO - ( x -11) I	 L - I	 1)- I	 I. - 1 (x - U
Y 1 D- I y
n
Y?	 di	 (5)
i=1
where
Y= L - 1 (X-U )	 ( 6)
Eq. (6) can be solved by forward substitution, i.e.
i-1
y i = (x i - u i ) -	 Iil y 	 (7)
j=1
where 1, = ( i i 1• ) . We then can use the Winograd procedure on the
2
summation term in eq. (7). (The staiJard method requires n - n
2
multiplies and (n- 2) 2( n - 1) adds to evaluate this term). We note
r
^ 	 l	 I	 I	 t	 l
8
that for evaluating this expression for more than one	 x , we can
use a variation of the algorithm for full matrix products.
Taking special note of the -,tructi : ,e of	 L,	 we then	 use
ri
i	 = L 'i, 2s- I	 t i, 2s
=1 i=:3, 4
	 ...	 ,	 n	 (8)
r.
^i	 = J y 2s-I y 2 s
s=1
i-	 I
with	 r i
	
=	 2 J Ibis	 is equivalent to
r.
- 3	= y y 2 j_2,;3,... 2	 (9)
In 2j	 "2j- I
TI2j+ I
	
^2j- 1
	
+ -`'.'j- 1
	 v 2j	 j=2, a, ... , ( n -
-	 2
We there have. with
i-I
a i	 L	 ^ij yi
j=l
a 2	 A21 yl
9ri
 
2j -I + y2j) (t
	
2j + y 2j-I ) -
j=1
R 	 - Ti	
if i odd
r	 r.
( t i ^ 2j-I
	
+ y2j) (ti, 2j + y2j-I)
j
F i	 ^i +	 t i, i- 1 y i- 1	 if	 i	 even
i=3 0 4, . .. , n	 (10)
Since ir depends only on I., it can be used in evaluating eq. (7)
for various x. "Fable 2 shows the number of operations necessary
for computing each of the various pars of the quadratic form. Note
that the Winograd algorithm is actually slower when Evaltiating this
expression for just one x, but for more than one x, the pre-
computed values of rr may be used. Table 3 shows the ratio of
times for the two methods for the case of 
"F in /T a = 2.7 for
computing ill of the quadratic terms for various values of n.
Note that a net savings occurs only in the lower right region.
Also shown in 'Fable 4 is the minimum ratio of T tn / T a for the
Winograd precedure to be faster for computing in quadratic terms.
Often in remote sensing calculations (e.g. maximum likelihood
classification), many probability density functions must be evaluated
over the same se ,., of data vectors. In this case, one may pre-
compute both the ^ 's and the I 's . Table 5 shows the asymp-
totic ( g precomputed and I computed only for the first class)
time savings to he exhccted for T  /7' a = 2.7 for k classes
for various values of n . Note that a net savings results for n > 4.
Ny
c0
H
N
b N
a O ^a^
E
"' o
U	 4"
C
CT —
w v
r- GIN
O1~
U N
^' CC
O
cz
Lo	 iIN
G w0 II
s. U	 R'
cz	
v
E U	 lr
v
z 3
G
bU
w
c0
%4
O
w
v
a
cl7
U
O
1
10
v
v
.ia
.r.
a-+
t-
uU
cU
t,
.O
.O
C
Ir
O CIN+
LL
NCI' rI , O C NCIN
cr
N
cr CT
+ d + N O O
cV
Q N Q
a a
m I N
N N+(^ CIN C C C
cn Id"
s~
Q"
M
^ ^^ C G
Q' iT
C-4
N
^
...
cr
c
R. +
N Cld a
E m N
•.r
cz
w ^
i00 N 00 Q0
N
O ? 00 OC 00
x r^ Q• Q, o•
N O O^ oG oo aoN
[- ^r
Q`
0
0
0
tTC
C- Q In N
OC C T Q` O^
Q+ N o0 to
Ln co CA CN
r0 Q^ O, C\
to N
N 0\ O% (31 Ol
00 to cry
N 00 00 00
O --^ O O, a, 0,
cn M :N r-+
x ^ C O O O
M ^:: V7 'r!' ^'f'
^O •-+ O O O O
N 00 I- t- 1-
rl" --^ O O O O
^^
O 0  OOt!7 O
11
w
O CJ
Vr
N
y L
V L1 ^
cz -r]
L ET7
W
RS	
^.L1 'C7
cz
h E `n o
O	 LrV
O Q1 .^
w ^
G
G
^ C v
O ._
V
O ^ ;
c
Q;
1 
. ,N.
`V
N
N
M
00
to 1 0
cn
w
In
N
N
I^
N
C14
CT
CA
.-.
[^
.--
N
'10
r i
all
••r
N
CV . r .•-i
N
r-
".4
N
 CT
•--^
C70
.--^
00
.-+
T
r-
N
en
cn
N
r^
cV N
/000V
.-+
N
I
p
C- 3
to
CV
M
N
.-.
N
r+
N
p
1
^i	 In
-3
M
1-
uo
11
%C
N
QO
-V
N
d'
N
7
N
•-+
<r M
00
p
M
Cn
CT
N
CT
QO
N
I	 ^0 to
In
O
I:
a0
M
r
C13
ti^ to
C14
to r-+ O
C14 i
VE^n
p
.-.
^n
c`1
O 0
O
1•r
u
w
O
P '^
0. wO
'C1
cti
cu
O
N 5 ^'
cd
co cQ'
E
E
-" bD
1~
.1
O cz
N U
M Cz
U
E G
2
^	 I	 I
13
^p cn .-.
^'T' to to 1n .n
N 0a 00 Cw) co
N T r- No
N to to LnN 00 00 00 GQ
clt 00 00 00 00
rJ, to M
n r,
r- oo 00 00 00
00 00 00 00
r 00 on 00 00
N h to m
^r O ON T O1
Q. 00 Cyr 00
N O O
.--+ T T T O^
^. co m N N
T T T T
T
00 %0 to ^r
T Q T ON
.c T 00 00 nON T T T
N N ^--^
IT O O O O
N Il — — 1
-4 ^^ I
ul
cc
.-r
U
^r
w
N^Octi
Z. E lc C
o -o z
4"co	 C
>, 1-4 0
 0Gov.
0r=^,^
G,^	 N
^ U	 1.4	 ?	
II
N
$4
^_ U O I-.
FG a U^
^v czcC U	 C
0 cd F-cz
Cr O `N
Q)
w l o
Q.
E cn
0 t
U
L ~
	 b^
O u^
t`^' O O
ca
U L+
O
E
cn
d
^ to O O
I	 1	 1	 ,^	 1
14
IV, Conclusions
'I'll( , 1A'inograd matrix multiplication sche , ,w can produce a time
savings in various computations using large order matrices in the
digital processing of remotely sensed data. We have shown ho w a
modification of this procedure may be applied to the calculation of
Gaussian probability density functions, and indicated 1ww this may be
extended to other computations. For large dimensions, or a large
numher of points, there can he some time savings, but the user
{	 should determine the expected time difference wing the value of
T m /T a of the computer to be used. A further study should be
undertaken to investigate the effects of decreased nurnerical stability
of this algorithm in v,.rious applications.
F
i
► 	 1 ,	 r
15
APPENDIX
SUNR7UTINc wNOMUL(A.P U•IA9R.SotB•C•IC• TMP I
C
C	 THIS 4OUT I NF USES TH I
 w INCG(?AD PPnC l7'DU Q1 T ]' MUL T I Q L Y Toll MA to I CF 5C	 t.F. C=A+ .j wHFPF-
 A 1S P * O: y . O + S: E C• P * S•
C	 IA IS TH= LFNGTH n F THE COLS OF A AS DI M ENSICN Fn IN TH= MAIN PROG•
C	 IN ^ IC AP P SIMILAR QUANTITIFS FOR D to C
C	 C MUS T NI T CrCUPY THE SAME STn2AGF AS A 14 H
TMP IS WIPKING STORA(, c7
 0= LFNGTH •GI7•S+P
t.
INTCGFG ) g Q.S•U• tJl.U2
PFAL A( IA•01.H( IBeSI * C.( IC.S)• T MP(I I
DOU31-^ P,a -CISMN SS.S1..2
L O GIC %L .IUD
I`TA=•P
L=Q/2
C
C	 COMPUTE THE X I f ^
C
00 10 1=1.P
IF (L.L T .2) 6 0 T l) 10
D r) 1'j U=2.L
U1=2*U
U'=U 1 —1
15	 SS=SS-A(I.U7)*A( I.UI )
10 TMP(I)=SS
C
C	 Cr'MPuTF TH c7 ETA I S
C
DO 20 J=1.S
:iS=-rn(1.J)+3(?.J)
IF (1_•L T.?) GO TD 20
10 2F- U=2 . L
U1=2+u
U2=U1-1
TMP( I F TA+J ) =SS
rDn=.FAL -, •
IF (2 + L.N P .0) IC)r)=•T(=UE.
r	 C
C	 C04PUTF THE C(I.J)•S
C
F Do 30 1 =1 . C
51=T'4 1-1 ( T )
S3=A( I.G)
n;l 30 J=I.
SS=Sl+TMP(I—;;TA+J)
IF(•IN rT * UDO1 G 	 TO .97
55=cS+S-4*H(0•J)
37	 Dr 3S U=1 .L
U1=2*u
U 2=11 l—1
35	 5	 5_;+(A(1.U2)+t,(U1.J))+(A(I.Ull+b(U?.J)!
10 C( I. J)=!^,S
FFTURN
FND
ORIGIN" PAGE V
OF POOR QUAUfyl
16
REFERENCES
[t ]	 Shmuel Winograd, - n the Number of Multiplications Necessary
to Compute Certain Functions," Comm. on Pure and Applied
Math. , 23, 1970, pp. 165-179.
[2] W. Miller, "Computational Complexity and Numerical Stability,"
SIAM J. Comput., 4, June, 1975, pp. 97-107.
[3] 1). L. Van Rooy, M. S. Lynn, any! C. 11. Snyder, "The Use
of the Modified Cholesky Decomposition in Divergence and
Classification Calculations, " ('.inf. Proc. Machine Processing;
of Remotely Sensed Data, LA AS, Rirduc, Oct. 16-18, 1973, 
pp. 313-35-313-47.


The use of the Winograd matrix multiplication algorithm in digital multispectral processing

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server