- Email: [email protected]
Fully-parallel relaxation algebraic operations for optical computers
Volume 43, number 4
OPTICS COMMUNICATIONS
15 October 1982
FULLY-PARALLEL RELAXATION ALGEBRAIC OPERATIONS FOR OPTICAL COMPUTERS ¢~
Wai K. CHENG * and H.J. CAULFIELD Aerodyne Research, Inc., The Research Center at Manning Park, Billerica, MA. 01821, USA
Received 22 June 1982
A new family of linear algebra algorithms for use with optical vector-times-matrixmultipliersis discussed. Rather than using iterative algorithms from the digital domain, we have used the equivalent differential analog implementation. The resuiting system relaxes continuously to the answer at a speed independent of the size of the matrix. These concepts are illustrated by a method for the solution of linear algebraic equations.
1. Introduction
Optical operations performable by properlyconnected optical matrix-times-vector operators are numerous [1-3]. These operators perform the operation for the N X N matrix A,
y2 Jill N
=A
x1 x2
(1)
x l N
or, equivalently, y =Ax
(2)
in parallel. That is all the N components o f y are calculated simultaneously. Typical use of these operators involves iterative operations [2,3 ] wherein (1) an input vectorx (0) is chosen, (2) y(0) = Ax(O) is calculated, (3) some operation is performed o n y (0) to calculate a new vector x (1), ( 4 ) y (1) = Ax (1) is calculated, and the process is repeated untily (k+l) andy (k) become essentially equal. These iterative algorithms are borrowed from digital computer domain. * Supported by Contracts N00014-82-C-0302 and F196288242-0068. * Consultant, also assistant professor of mechanicalengineering, Massachusetts Institute of Technology. 0 030-4018/82/0000-0000/$ 02.75 © 1982 North-Holland
Implementation of these algorithms often use a hybrid combination of the optical matrix multiplier and a digital computer which controls the data path and steps through the iteration. Furthermore, most applications of those algorithms have not been fully parallel. Sequential treatments of the parallel data from the optical operators are processed by the digital computer in the loop. Accordingly, our goals have been twofold: (1) design fully parallel algorithms and hardware, (2) design analog feedback (relaxation) implementations rather than the iterative procedures. These fully-parallel, relaxation algorithms should converge to the proper answer at a speed which is independent of the vector length N. What governs the convergence speed is (1) the inherent hardware speed and (2) the problem details. In what follows, we will explain one algorithm in detail. Others are easily derived by analogy and will be presented in subsequent publications. Our purpose is to show a principle as well as to illustrate the principle with a specific example.
2. Solution to simultaneous equations
In this example we shall consider the problem of solving for the vector x from the equation A x =y ,
(3) 251
Volume 43, number 4
OPTICS COMMUNICATIONS (Y -
15 October 1982
A.x) 1
N [ y . ei ~
-I
x ( t - + ° ° ) = i ~ 1 ~--~-i ) e i = A
(y
-
If the real parts of Xi are all negative, the signs of the difference amplifiers D 1 and D 2 may be reversed and the correct asymptotic solution may still be obtained. If there are one or more of the real part of the eigenvalues of opposite signs, however, the solutions always diverge. This divergence will cause saturation of the amplifiers. To guarantee stability of the system, the feedback system of fig. l is modified to implement the following transfer function:
Ax) 2
Fig. 1. A conditionally stable matrix equation solver for a
2 × 2 matrix. wherey is a given realN X 1 vector. The N × N real matrix A is assumed to be of full rank and may have eigenvalues which are in general complex and may not be distinct. For purpose of illustration, let A be 2 X 2. The implementation can easily be generalized to the case of a n n X N matrix. Consider first the circuit implementation as shown in fig. 1. The vectory is represented by the voltage sources Y 1 and Y2" The construction of the circuit guarantees that if a stable equilibrium is attended, the resulting voltages x 1 and x 2 will be the solution. Nevertheless, it will be shown that the stability of the circuit depends on the eigenvalues of the matrix A. The stability of the circuit may be analyzed by looking at its asymptotic behavior for an arbitrary starting condition. For this case, it suffices to assume distinct eigenvalues for A. The circuit is described by the differential equation
(d / dr) x = -(1/r)(Ax - y),
(8)
where D =
"..
.
(9)
XN
The solution to (7) for an arbitrary initial condition x 0 is then
x=xT~exp(-Xlt)
0 ]
"..
0
Sx 0
exp(-X n t)
I - e x p ( - X l t)
+ ST
I
Xl
0
1
SATy(IO)
• 1 - exp(-Xnt!l
0
v.
4
(1 - exp(-~tit))J ei,
(5)
where the e i are the eigenvectors with )ki the corresponding eigenvalues. If all the eigenvalues Xi are real and positive, the solution of (5) indeed converges to 252
S A T A S T = D,
-j
The positive definiteness of ?ti guarantees stability of the system and the asymptotic solution o f x is
i=1 L
Xi
(7)
Here the superscript T indicates the matrix transpose. The corresponding circuitry is shown in fig. 2, Since ATA is a real and symmetric matrix all its eigenvalues h i are real and positive, and there exists a real orthogonal matrix S such that
i-
X(t) = ~ Ix 0 • e i exp(-Xit) +Y" e i
(d/dt)x = - ( 1 / r ) A T ( A x - y).
(4)
where r is the time response of the integrators I 1 and 12. We may set r equal to 1 by measuring t in units of r. The solution to (4), for an initial condition x0, is N
(6)
Y.
x(t ~°°)
ST
o1
SATy.
•
(11)
1
qj
From (8), postmultiplying by (AST) -1 , we obtain
Volume 43, number 4
OPTICS COMMUNICATIONS [AT(Ax- Y)"] l [AT(Ax - ~)] 2
l'
out
AT
in
15 October 1982
I o,I\
/+
\ D2
11 ) y1 |_
Fig. 2. A stable matrix equation solver. (12)
SA T =DSA -1 .
Combining (11) and (12) we obtain the desired result (13)
x =AMy.
The convergence o f (8) from an arbitrary initial condition x 0 therefore depends on the minimum value Xm o f the set of h I and goes as e x p ( - k m t ) regardless of the size of the matrix. This is a result of the fully parallel nature of the implementation. The above argument can be applied to a complex matrix A and a complex vectory. Then complex number multipliers have to be used and the transpose matrices A T, S T will be replaced by the hermitian adjoints A t and S t . The matrix S will then be unitary. When the matrix A is not of full rank, the matrix A T A will be positive semi-definite. The algorithm will relax to a vector x which satisfies (3). Since x is not unique, the particular asymptotic value for x will be dependent on the initial condition x 0. A reader of our manuscript wondered why the magnitude of ~ could be outside the unit circle for this continuoUs relaxation method, but is restricted to within the unit circle for the discrete case. That question has an interesting answer. In the discrete case there is a single clock time during which all operations take place. In the continuous case dx/dt = --(1/r) A T A x
(14)
implies x (n+l) -- x (n) At
1 - -- X2x(n), r
(15)
where x (n+l) and x (n) are samples taken at times At apart and k 2 > 0. Because, unlike the discrete case, we can not equate At with r, the convergence requirement is
1 -&--tx21r <1
(16)
not, as before, I1 -X21 < 1.
(17)
But r is fixed and At is arbitrarily small. Thus there is no maximum allowable value of IX I. The same reader also pointed out that the use of A T A to assure a positive eigenvalue had been used earlier this year by Goodman and Song [4]. We are grateful to Optics Communications and its excellent reader for asking such a helpful question and pointing out such an important omitted reference.
3. Discussion The relaxation algorithm we have described is clearly the differential version of the iterative algorithm [2,3]. Our "algorithm" for generating these relaxation algorithms is to drive the circuit in such a way as to embody the definition of the solution by driving the inputs using the error signals. In this case we drive x using A x - y . There is no definite number of iterations. Rather A x - y converges smoothly (not in discrete jumps) to zero. This approach appears to
253
Volume 43, number 4
OPTICS COMMUNICATIONS
take full advantage of the continuous nature of the optical processor and to require that continuous nature to assure convergence for eigenvalues on or outside the unit circle.
References [1 ] J.W. Goodman, A.R. Dias and L.M. Woody, Optics Lett. 2 (1978) 1.
254
15 October 1982
[2] D. Psaltis, D. Casasent and M. Carlotto, Optics Lett. 4 (1979) 348. [3] H.J. Caulfield, D. Dvore, J.W. Goodman and W.H. Rhodes, Appl. Optics 20 (1981) 2263. [4] J.W. Goodman and M.S. Song, Appl. Optics 21 (1982) 502.
OPTICS COMMUNICATIONS
15 October 1982
FULLY-PARALLEL RELAXATION ALGEBRAIC OPERATIONS FOR OPTICAL COMPUTERS ¢~
Wai K. CHENG * and H.J. CAULFIELD Aerodyne Research, Inc., The Research Center at Manning Park, Billerica, MA. 01821, USA
Received 22 June 1982
A new family of linear algebra algorithms for use with optical vector-times-matrixmultipliersis discussed. Rather than using iterative algorithms from the digital domain, we have used the equivalent differential analog implementation. The resuiting system relaxes continuously to the answer at a speed independent of the size of the matrix. These concepts are illustrated by a method for the solution of linear algebraic equations.
1. Introduction
Optical operations performable by properlyconnected optical matrix-times-vector operators are numerous [1-3]. These operators perform the operation for the N X N matrix A,
y2 Jill N
=A
x1 x2
(1)
x l N
or, equivalently, y =Ax
(2)
in parallel. That is all the N components o f y are calculated simultaneously. Typical use of these operators involves iterative operations [2,3 ] wherein (1) an input vectorx (0) is chosen, (2) y(0) = Ax(O) is calculated, (3) some operation is performed o n y (0) to calculate a new vector x (1), ( 4 ) y (1) = Ax (1) is calculated, and the process is repeated untily (k+l) andy (k) become essentially equal. These iterative algorithms are borrowed from digital computer domain. * Supported by Contracts N00014-82-C-0302 and F196288242-0068. * Consultant, also assistant professor of mechanicalengineering, Massachusetts Institute of Technology. 0 030-4018/82/0000-0000/$ 02.75 © 1982 North-Holland
Implementation of these algorithms often use a hybrid combination of the optical matrix multiplier and a digital computer which controls the data path and steps through the iteration. Furthermore, most applications of those algorithms have not been fully parallel. Sequential treatments of the parallel data from the optical operators are processed by the digital computer in the loop. Accordingly, our goals have been twofold: (1) design fully parallel algorithms and hardware, (2) design analog feedback (relaxation) implementations rather than the iterative procedures. These fully-parallel, relaxation algorithms should converge to the proper answer at a speed which is independent of the vector length N. What governs the convergence speed is (1) the inherent hardware speed and (2) the problem details. In what follows, we will explain one algorithm in detail. Others are easily derived by analogy and will be presented in subsequent publications. Our purpose is to show a principle as well as to illustrate the principle with a specific example.
2. Solution to simultaneous equations
In this example we shall consider the problem of solving for the vector x from the equation A x =y ,
(3) 251
Volume 43, number 4
OPTICS COMMUNICATIONS (Y -
15 October 1982
A.x) 1
N [ y . ei ~
-I
x ( t - + ° ° ) = i ~ 1 ~--~-i ) e i = A
(y
-
If the real parts of Xi are all negative, the signs of the difference amplifiers D 1 and D 2 may be reversed and the correct asymptotic solution may still be obtained. If there are one or more of the real part of the eigenvalues of opposite signs, however, the solutions always diverge. This divergence will cause saturation of the amplifiers. To guarantee stability of the system, the feedback system of fig. l is modified to implement the following transfer function:
Ax) 2
Fig. 1. A conditionally stable matrix equation solver for a
2 × 2 matrix. wherey is a given realN X 1 vector. The N × N real matrix A is assumed to be of full rank and may have eigenvalues which are in general complex and may not be distinct. For purpose of illustration, let A be 2 X 2. The implementation can easily be generalized to the case of a n n X N matrix. Consider first the circuit implementation as shown in fig. 1. The vectory is represented by the voltage sources Y 1 and Y2" The construction of the circuit guarantees that if a stable equilibrium is attended, the resulting voltages x 1 and x 2 will be the solution. Nevertheless, it will be shown that the stability of the circuit depends on the eigenvalues of the matrix A. The stability of the circuit may be analyzed by looking at its asymptotic behavior for an arbitrary starting condition. For this case, it suffices to assume distinct eigenvalues for A. The circuit is described by the differential equation
(d / dr) x = -(1/r)(Ax - y),
(8)
where D =
"..
.
(9)
XN
The solution to (7) for an arbitrary initial condition x 0 is then
x=xT~exp(-Xlt)
0 ]
"..
0
Sx 0
exp(-X n t)
I - e x p ( - X l t)
+ ST
I
Xl
0
1
SATy(IO)
• 1 - exp(-Xnt!l
0
v.
4
(1 - exp(-~tit))J ei,
(5)
where the e i are the eigenvectors with )ki the corresponding eigenvalues. If all the eigenvalues Xi are real and positive, the solution of (5) indeed converges to 252
S A T A S T = D,
-j
The positive definiteness of ?ti guarantees stability of the system and the asymptotic solution o f x is
i=1 L
Xi
(7)
Here the superscript T indicates the matrix transpose. The corresponding circuitry is shown in fig. 2, Since ATA is a real and symmetric matrix all its eigenvalues h i are real and positive, and there exists a real orthogonal matrix S such that
i-
X(t) = ~ Ix 0 • e i exp(-Xit) +Y" e i
(d/dt)x = - ( 1 / r ) A T ( A x - y).
(4)
where r is the time response of the integrators I 1 and 12. We may set r equal to 1 by measuring t in units of r. The solution to (4), for an initial condition x0, is N
(6)
Y.
x(t ~°°)
ST
o1
SATy.
•
(11)
1
qj
From (8), postmultiplying by (AST) -1 , we obtain
Volume 43, number 4
OPTICS COMMUNICATIONS [AT(Ax- Y)"] l [AT(Ax - ~)] 2
l'
out
AT
in
15 October 1982
I o,I\
/+
\ D2
11 ) y1 |_
Fig. 2. A stable matrix equation solver. (12)
SA T =DSA -1 .
Combining (11) and (12) we obtain the desired result (13)
x =AMy.
The convergence o f (8) from an arbitrary initial condition x 0 therefore depends on the minimum value Xm o f the set of h I and goes as e x p ( - k m t ) regardless of the size of the matrix. This is a result of the fully parallel nature of the implementation. The above argument can be applied to a complex matrix A and a complex vectory. Then complex number multipliers have to be used and the transpose matrices A T, S T will be replaced by the hermitian adjoints A t and S t . The matrix S will then be unitary. When the matrix A is not of full rank, the matrix A T A will be positive semi-definite. The algorithm will relax to a vector x which satisfies (3). Since x is not unique, the particular asymptotic value for x will be dependent on the initial condition x 0. A reader of our manuscript wondered why the magnitude of ~ could be outside the unit circle for this continuoUs relaxation method, but is restricted to within the unit circle for the discrete case. That question has an interesting answer. In the discrete case there is a single clock time during which all operations take place. In the continuous case dx/dt = --(1/r) A T A x
(14)
implies x (n+l) -- x (n) At
1 - -- X2x(n), r
(15)
where x (n+l) and x (n) are samples taken at times At apart and k 2 > 0. Because, unlike the discrete case, we can not equate At with r, the convergence requirement is
1 -&--tx21r <1
(16)
not, as before, I1 -X21 < 1.
(17)
But r is fixed and At is arbitrarily small. Thus there is no maximum allowable value of IX I. The same reader also pointed out that the use of A T A to assure a positive eigenvalue had been used earlier this year by Goodman and Song [4]. We are grateful to Optics Communications and its excellent reader for asking such a helpful question and pointing out such an important omitted reference.
3. Discussion The relaxation algorithm we have described is clearly the differential version of the iterative algorithm [2,3]. Our "algorithm" for generating these relaxation algorithms is to drive the circuit in such a way as to embody the definition of the solution by driving the inputs using the error signals. In this case we drive x using A x - y . There is no definite number of iterations. Rather A x - y converges smoothly (not in discrete jumps) to zero. This approach appears to
253
Volume 43, number 4
OPTICS COMMUNICATIONS
take full advantage of the continuous nature of the optical processor and to require that continuous nature to assure convergence for eigenvalues on or outside the unit circle.
References [1 ] J.W. Goodman, A.R. Dias and L.M. Woody, Optics Lett. 2 (1978) 1.
254
15 October 1982
[2] D. Psaltis, D. Casasent and M. Carlotto, Optics Lett. 4 (1979) 348. [3] H.J. Caulfield, D. Dvore, J.W. Goodman and W.H. Rhodes, Appl. Optics 20 (1981) 2263. [4] J.W. Goodman and M.S. Song, Appl. Optics 21 (1982) 502.