High Precision Computing With Charge Domain Devices and a Pseudo-Spectral Method Therefor

1. Origin of the Invention

The invention described herein was made in the performance of work under a NASA contract, and is subject to the provisions of Public Law 96-517 (35 USC 202) in which the contractor has elected to retain title.

2. Technical Field

The invention relates to charge coupled device (CCD) and charge injection device (CID) hardware applied to higher precision parallel arithmetic processing devices particularly adapted to perform large numbers of multiply-accumulate operations with massive parallelism and to methods for solving partial differential equations therewith.

3. Background Art

"Grand Challenges" have been defined as fundamental problems in science or engineering, with broad economic and scientific impact, that could be advanced by applying high performance computing resources, while high speed digital computers with some level of parallelism are enabling the consideration of an ever growing number of practical applications, many very large-scale problems await the development of massively parallel hardware. To carry out the demanding computations involved in grand challenges, it is generally accepted that one needs to pursue a hybrid approach involving both novel algorithms development and revolutionary chip technology.

Many algorithms required for scientific modeling make frequent use of a few well defined, often functionally simple, but computationally very intensive data processing operations. Those operations generally impose a heavy burden on the computational power of a conventional general-purpose computer, and run much more efficiently on special-purpose processors that are specifically tuned to address a single intensive computation task only. A typical example among the important classes of demanding computations are vector and matrix operations such as multiplication of vectors and matrices, solving linear equations, matrix inversion, eigenvalue and eigenvector search, etc. Most of the computationally more complex vector and matrix operations can be reformulated in terms of basic matrix-vector and matrix-matrix multiplications. From a neutral network perspective, the product of the synaptic matrix by the vector of neuron potentials is another good example.

An innovative hybrid, analog-digital charge-domain technology, for the massively parallel VLSI implementation of certain large scale matrix-vector operations, has recently been developed, as disclosed in U.S. Pat. No. 5,054,040. It employs arrays of Charge Coupled/Charge Injection Device (CCD/CID) cells holding an analog matrix of charge, which process digital vectors in parallel by means of binary, non-destructive charge transfer operations. FIG. 1 shows a simplified schematic of the CCD/CID array 100. Each cell 110 in the array 100 connects to an input column line 120 and an output row line 130 by means of a column gate 140 and a row gate 150. The gates 140, 150 hold a charge packet 160 in the silicon substrate underneath them that represents an analog matrix element. The matrix charge packets 160 are initially stored under the column gates 140. In the basic matrix-vector multiplication (MVM) mode of operation, for binary input vectors, the matrix charge packets 160 are transferred from under the column gate 140 toward the row gates 150 only if the input bit of the column indicates a binary `one`. The charge transferred under the row gates 150 is summed capacitively on each output row line 130, yielding an analog output vector which is the product of the binary input vector with the analog charge matrix. By virtue of the CCD/CID device physics, the charge sensing at the row output lines 130 is of a non-destructive nature, and each matrix charge packet 160 is restored to its original state simply by pushing the charge back under its column gate 140.

FIGS. 2A, 2B, 2C, 2D and 2E are an illustration of the binary-analog MVM computation cycle for a single row of the CCD/CID array. In FIG. 2A, the matrix charge packet 160 sits under the column gate 140. In FIG. 2B, the row line 130 is reset to a reference voltage. In FIG. 2C, if the column line 120 receives a logic one input bit, the charge packet 160 is transferred underneath the row gate 150. In FIG. 2D, the transferred charge packet 160 is sensed capacitively by a change in voltage on the output row line 130. In FIG. 2E the charge packet 160 is returned under the column gate 140 in preparation for the next cycle. A bit-serial digital-analog MVM can be obtained from a sequence of binary-analog MVM operations, by feeding in successive vector input bits sequentially, and adding the corresponding output contributions after scaling them with the appropriate powers of two. A simple parallel array of divide-by-two circuits at the output accomplishes this task. Further extensions of the basic MVM scheme of FIG. 1 support full digital outputs by parallel A/D conversion at the outputs, and four-quadrant operation by differential circuit techniques.

FIGS. 3A and 3B illustrate how the matrix charge packets are loaded into the array. In FIG. 3A, appropriate voltages are applied to each gate 170 in each cell 110 of the CCD/CID array 100 so as to configure each cell 110 as a standard 4-phase CCD analog shift register to load all of the cells 110 sequentially. In FIG. 3B, the same gates 170 are used for row and column charge transfer operations as described above with reference to FIG. 2.

The particular choice for this unusual charge-domain technology resulted from several considerations, not limited to issues of speed and parallelism. In comparison to other, more common parallel high-speed technology environments (digital CMOS, etc.), the distinct virtues of the foregoing charge-domain technology for large-scale special-purpose MVM operations are the following: .smallcircle. Very High Density: The compactness of the CCD/CID cell allows the integration of up to 10.sup.5 on a 1 cm.sup.2 die (in a standard 2 .mu.m CMOS technology), providing single-chip 100 GigaOPS computation power. .smallcircle. Very Low Power Consumption: The charge stored in the matrix is conserved along the computation process because of the non-destructive nature of the CCD/CID operation. Hence, the entire power consumption is localized at the interface of the array, for clocking, I/O and matrix refresh purposes. This enables the processor to operate at power levels in the mW/TeraOPS range. .smallcircle. Scalability: The scalable architecture of the CCD/CID array allows the interfacing of many individual processors in parallel, combining together to form effective processing units of higher dimensionality; still operating at nominal speed. .smallcircle. I/O Flexibility: Although an analog representation is used inside the array to obtain fast parallel computation, the architecture of the processor provides the flexibility of a full digital interface, eliminating the bandwidth and noise problems typical for analog interfacing. .smallcircle. Programming Flexibility: The architecture allows for either optical (parallel, sustained) or electronic (semi-parallel, periodical) loading of the charge matrix. The latter method, described above with reference to FIG. 3, requires interrupts, of duration usually much shorter than the time interval in computation mode, for which the stored charge matrix remains valid before a matrix refresh is needed.

Preliminary results on a 128.times.128 working prototype, implemented on a single 4 mm.times.6 mm die in 2 .mu.CCD-CMOS technology, indicate a performance of approximately 10.sup.10 8-bit multiply-accumulate operations per second. Processors with 1024.times.1024 cells will be realizable in the near future.

An innovative approach to parallel algorithms design remains a key enabling factor for high performance computing. In contrast to conventional approaches, one must develop computational schemes which exploit, from the onset, the concept of massive parallelism. The CCD/CID MVM approach discussed above represents a promising hardware technology for massively parallel computation, which offers both opportunities and challenges in the design of parallel algorithms. This new technology defies some of the most basic assumptions in the analysis of computational complexity of parallel algorithms. To clarify this, a short discussion is now given. In what follows, N.times.1 vectors and N.times.N matrices are considered. Also, m and a stand for multiplication and addition, respectively. One of the most fundamental results regarding the complexity of parallel computation involves the summation of N scalars.

The summation of N scalars can be performed in O(log N).alpha. operations with O(N) processors. The complexity of some other matrix-vector operations follows from the above theorem as:

Corollary 1.1. A vector-dot product can be performed in O(log N).alpha.+O(1)m operations with O(N) processors.

Corollary 1.2. A matrix-vector multiplication can be performed in O(log N).alpha.+O(1)m operations with O(N.sup.2) processors.

Corollary 1.3. A matrix-matrix multiplication can be performed in O(log N).alpha.+O(1)m operations with O(N.sup.3) processors.

The complexity of many other computational problems, including Fourier transforms and linear systems, can also be derived based on the above theorem and its corollaries. Currently, these results provide the framework for the design of efficient parallel algorithms. It should be mentioned that Corollaries 1.2 and 1.3 are more of a theoretical importance than a practical one, since the implementation of parallel algorithms achieving the above time-lower bounds requires rather complex two-dimensional processor interconnections. Interestingly, the result of Theorem 1 has been assumed to be independent of hardware technology. This reflects a basic feature of digital computers, which do not allow the simultaneous summation of N numbers.

The CCD/CID technology of FIGS. 1-3 defies this most basic assumption, since the summation of N numbers is just the summation of N charges, which can be performed simultaneously. This is a fundamental change, which leads to a new set of results regarding the complexity of parallel computation using CCD/CID architectures.

The summation of N scalars can be performed in O(1).alpha. operations with O(N) "processors" (CCD/CID cells). The complexity of matrix-vector operations follows from the above theorem. In particular:

Corollary 2.1. A vector-dot product can be performed in O(1).alpha.+O(1)m operations with O(N) "processors" or CCD/CID cells.

Corollary 2.2. A matrix-vector multiplication can be performed in O(1.alpha.+O(1)m operations with O(N.sup.2) "processors" or CCD/CID cells.

Corollary 2.3. A matrix-matrix multiplication can be performed in O(1).alpha.+O(1)m operations with O(N.sup.3) "processors" or CCD/CID cells.

The complexity of these operations is independent of the problem size, assuming that the chip size matches the problem. The above results are not only significant from a theoretical point of view, since they improve the known time-lower bounds in these computations, but are also of practical importance. They show that the design and analysis of parallel algorithms based on CCD/CID architectures is drastically difference from that for conventional parallel computers. The CCD/CID chip is currently considered for performing Matrix-Vector Multiplication (MVM) for which it achieves the computation time given in Corollary 2.2. However, in order to output the results, k clock cycles are needed, where k is the required precision, in number of bits. The CCD/CID chip is most efficient for MVM computations wherein the matrix is known beforehand, as well as for series of MVMs performed with the same matrix. Despite this seemingly narrow function, the CCD/CID processor can be used for many applications, leading to new results. In the following, a few examples are given:

a. The Discrete Fourier Transform (DFT). The DFT can be presented as a MVM for which the time-lower bound can be derived from Corollary 2.2. The DFT represents an application for which the coefficient matrix is known beforehand. However, for both serial and conventional parallel computation, the FFT is always preferred due to its greater efficiency. In particular, for parallel computation, the asymptotic time lower bound of O(log.sub.2 N) can be achieved by using O(N) "processors". In contrast, for the CCD/CID chip the DFT is more efficient for large problems, and can be performed in O(k) steps with O(N.sup.2) CCD/CID cells. The computational complexity of O(k) not only improves the previously assumed time-lower-bound, but is also independent of the transformed size.

b. Linear System Solution. For many applications, the solution of a linear system, Ax=b, is required, wherein the coefficient matrix A is known a priori. This is, for example, the case for the solution of partial differential equations (PDE) using a finite difference method, e.g., parabolic PDE for which A is a tridiagonal matrix, or Poisson equation, for which A is block tridiagonal. For both cases, given the fact that A is positive definite, its inverse A.sup.-1 can also be computed a priori. The linear systems solution can then be reduced to a simple MVM. However, this is not efficient for either serial or conventional parallel computation since, unlike A, the matrix A.sup.-1 is dense. There are more efficient algorithms for both serial and parallel computation which exploit the sparse structure of A. For the CCD/CID chip, on the other hand, the unconventional and seemingly inefficient method is more suitable, since it involves the above linear system in O(k) steps with O(N.sup.2) CCD/CID cells.

c. Eigenvector Search. Given a matrix A having N linearly independent eigenvectors and associated eigenvalues, one often needs to find the eigenvector v.sup.(1) corresponding to the dominant eigenvalue .lambda..sub.1. One available technique is the Power Method which, starting from an arbitrary vector u.sup.(o), performs successive MVMs u.sup.(n) =Au.sup.(n-1). With a suitable scaling, this sequence converges to cv.sup.(1), for some constant c.

It is believed that charge-domain processors not only represent a new hardware computing technology, allowing high performance solution for a large class of problems, but that they also affect many theoretical issues regarding the design and analysis of massively parallel algorithms. One area of particular importance for future advances is the need to extend the device's dynamic range (accuracy).

The advances in analog VLSI architectures discussed above provide a strong incentive to develop massively parallel algorithms for selected information processing problems. However, even though the hardware available to date is characterized by very high speed, its precision is typical of CCD/CID technology, i.e., of the order of 8 bits. Since each matrix element in FIG. 1 is represented by a charge packet, its accuracy is limited to the resolution of a CCD charge packet. A typical CCD charge packet is on the order of a million electrons, and this size limits its resolution to under 8 bits. Furthermore, voltage resolution (i.e., charge sensing) at the end of each row line is also of the order of 8 to 10 bits. Thus, both architectures and algorithms are needed, which not only are able to fully utilize the hardware's throughput, but also, ultimately, provide highly accurate results. Accordingly, it is an object of the present invention to extend the accuracy or resolution of CCD/CID MVM processors far beyond that currently available. It is a related object of the invention to do so without incurring a proportionate penalty in speed.

The present invention enhances the bit resolution of a CCD/CID MVM processor by storing each bit of each matrix element as a separate CCD charge packet. The bits of each input vector are separately multiplied by each bit of each matrix element in massive parallelism and the resulting products are combined appropriately to synthesize the correct product. In one embodiment, the CCD/CID MVM array is a single planar chip in which each matrix element occupies a single column of b bits, b being the bit resolution, there being N rows and N columns of such single columns in the array. In another embodiment, the array constitutes a stack of b chips, each chip being a bit-plane and storing a particular significant bit of all elements of the matrix. In this second embodiment, an output chip is connected edge-wise to the bit-plane chips and performs the appropriate arithmetic combination steps. Such arrays are employed in a pseudo-spectral method of the invention, in which partial differential equations are solved by expressing each derivative analytically as matrices, and the state function is updated at each computation cycle by multiplying it by the matrices. The matrices are treated as synaptic arrays of a neural network and the state function vector elements are treated as neurons. Using such a neural architecture, both linear (e.g., heat) and non-linear (e.g., soliton) partial differential equations are solved. Moving target detection is performed by driving the soliton equation with a vector of detector outputs. The neural architecture consists of two synaptic arrays corresponding to the two differential terms of the soliton equation and an adder connected to the output thereof and to the output of the detector array to drive the soliton equation. A resonance phenomenon between the moving waves of the soliton equation and moving targets sensed by the detector array is observed which enhances the amplitude of the state function neuron elements at locations corresponding to target locations.

In a preferred embodiment, the MVM processor of the invention includes an array of N rows and M columns of CCD matrix cell groups corresponding to a matrix of N rows and M columns of matrix elements, each of the matrix elements representable with b binary bits of precision, each of the matrix cell groups including a column of b CCD cells storing b CCD charge packets representing the b binary bits of the corresponding matrix element, the amount of charge in each packet corresponding to one of two predetermined amounts of charge. Each of the CCD cells includes a holding site and a charge sensing site, each charge packet initially residing at the respective holding site. The MVM processor further includes a device for sensing, for each row, an analog signal corresponding to a total amount of charge residing under all charge sensing sites of the CCD cells in the row, an array of c rows and M columns CCD vector cells corresponding to a vector of M elements representable with c binary bits of precision, each one of the M columns of CCD vector cells storing a plurality of c charge packets representing the c binary bits of the corresponding vector element, the amount of charge in each packet corresponding to one of two predetermined amounts of charge. A multiplying device operative for each one of the c rows of the CCD vector cells temporarily transfers to the charge sensing site the charge packet in each one of the M columns of matrix cells for which the charge packet in the corresponding one of the M columns and the one row of the CCD vector cells has an amount of charge corresponding to a predetermined binary value.

The preferred embodiment further includes an arithmetic processor operative in synchronism with the multiplying a device including a device for receiving, for each row, the sensed signal, whereby to receive N.times.b signals in each one of c operations of the multiplying a device, a device for converting each of the signals to a corresponding byte of output binary bits, and a device for combining the output binary bits of all the signals in accordance with appropriate powers of two to generate bits representing an N-element vector corresponding to the product of the vector and the matrix.

In one embodiment, the array of matrix CCD cells is distributed among a plurality of b integrated circuits containing sub-arrays of the M columns and N rows of the matrix CCD cells, each of the sub-arrays corresponding to a bit-plane of matrix cells representing bits of the same power of two for all of the matrix elements. A backplane integrated circuit connected edgewise to all of the b integrated circuits includes a device for associating respective rows of the vector CCD elements with respective rows of the matrix CCD elements, whereby the multiplying device operates on all the rows of the vector CCD elements in parallel.

The invention further includes a neural architecture for solving a partial differential equation in a state vector consisting of at least one term in a spatial partial derivative of the state vector and a term in a partial time derivative of the state vector. The neural architecture includes a device for storing an array of matrix elements of a matrix for each partial derivative term of the equation, the matrix relating the partial derivative to the state vector, a device for multiplying all rows of the matrix elements by corresponding elements of the state vector simultaneously to produce terms of a product vector, and a device for combining the product vector with a previous iteration the state vector to product a next iteration of the state vector. The equation can include plural spatial partial derivative terms wherein the neural architecture includes plural arrays of matrix elements of matrices corresponding to the partial spatial derivative terms and a device for combining the product vectors. The partial spatial derivative terms are definable over a spatial grid of points and the matrix elements are obtained in a pseudo-spectral analytical expression relating each matrix element to the distance between a corresponding pair of rid points.

In one embodiment of the neural architecture, moving targets are detected in the midst of noise in images sensed by an array of sensors, by the addition to the sum of product vectors of a vector derived from the outputs of the sensors, to produce a next iteration of the state vector. Preferably, the stored matrices implement a soliton equation having first and third order partial spatial derivatives of the state vector, so that the neural architecture includes a device for storing respective arrays corresponding to the first and third order derivatives. This produces two product vectors corresponding to the first and third order derivatives, respectively. This embodiment further includes a device for multiplying the product vector of the first order derivative by the state vector.

FIG. 1 is a simplified diagram of a CCD/CID MVM processor array of the prior art.

FIGS. 2A, 2B, 2C, 2D and 2E illustrate a sequence of matrix-vector multiply operations in a unit cell of the array of FIG. 1.

FIGS. 3A and 3B illustrate, respectively, electronic loading of the matrix elements and arithmetic operations in a unit cell of the array of FIG. 1.

FIG. 4 is a plan view of a CCD/CID MVM processor of the present invention.

FIG. 5 is a schematic diagram of a typical architecture for a higher precision arithmetic processor employed in combination with the CCD/CID processor of FIG. 4.

FIG. 6 is a diagram of a three-dimensional embodiment of the CCD/CID MVM processor of the present invention.

FIG. 7 is a schematic diagram of a general neural architecture for implementing the pseudo-spectral method of the present invention applied to linear partial differential equations.

FIG. 8 is a schematic diagram of a neural architecture for performing the pseudo-spectral method applied to a non-linear partial differential equation and for performing target detection therewith.

FIGS. 9A-9D include data plots obtained in a computer simulation of the neural architecture of FIG. 8 in which the sensors only sense noise, of which FIG. 9A is a Fourier spectrum of the sensor input, FIG. 9C is the network output and FIG. 9D is the filtered network output in which noise has been removed.

FIGS. 10A-10D include data plots in a computer simulation of the neural architecture of FIG. 8 in which the sensors sense a moving target with a per-pulse signal-to-noise ratio of -0 dB, of which FIG. 10A is a Fourier spectrum of the sensor input, FIG. 10B is a spatial plot of the sensor inputs, FIG. 10C is the network output and FIG. 10D is the filtered network output in which noise has been removed.

FIGS. 11A-11D include data plots obtained in a computer simulation of the neural architecture of FIG. 8 in which the sensors sense a moving target with a per-pulse signal-to-noise ratio of -10 dB, of which FIG. 11A is a Fourier spectrum of the sensor input, FIG. 11B is a spatial plot of the sensor inputs, FIG. 11C is the network output and FIG. 11 is the filtered network output in which noise has been removed.

In order to achieve high precision in a CCD/CID MVM processor, the present invention provides an advanced architecture which is illustrated schematically in FIGS. 4 and 5. In the following, the invention is described in its basic form. However, different architectures can be derived from this basic from which are not discussed here. A key element of the CCD/CID MVM array 200 of FIGS. 4 and 5 is the encoding in each CCD/CID processor cell 210 of one bit of the binary representation of each matrix element. As shown in FIG. 4, if a matrix A is to be specified with b bits of precision, each element A.sub.ij of the matrix occupies a single column 205 of b cells, namely cells corresponding to [b(i-1]+l,j], for l=1 , . . . b, labelled in FIG. 4 as A.sub.ij.sup.0 , . . . A.sub.ji.sup.b-1. Thus, the CCD/CID MVM processor array 200 of FIG. 4 is an array of N columns and N rows of matrix elements, each matrix element itself being a column 205 of b CCD/CID cells 210. There are, therefore, at total of N.times.N.times.b cells 210 in the array 200. A vector of N elements representable with a precision of b binary bits is stored in an array 230 of CCD/CID cells 235, the array 230 being organized in N columns 240 of b CCD/CID cells, each cell 235 in a given column 240 storing the corresponding one of the b binary bits of the corresponding vector element. Each CCD/CID cell 210 in FIG. 4 is of the type described above with reference to FIGS. 1-3 and is operated in the same manner with N column input lines (coupled to a successive row of N vector CCD/CID cells 235 in the vector CCD/CID array 230) and row output lines of the type illustrated in FIG. 1.

Matrix-vector multiplication will now be described with reference to the example of a conventional matrix-vector product obtained by adding the products of the elements in the vector and each row of the matrix. Of course, the present invention is not confined to a particular type of matrix vector product, and can be used to compute other types of matrix-vector products. One example of anther well-known type of matrix-vector product is that obtained by adding the products of a respective vector element multiplied by all elements in a respective column of the matrix.

Computation proceeds as follows. At clock cycle one, the matrix A, in its binary representation, is multiplied by the binary vector labelled u.sub.1.sup.0 , . . . u.sub.N.sup.0, which contains the least significant bits of u.sub.1 , . . . u.sub.N (i.e., by the top row of charge packets in the array 230 of vector CCD/CID cells 235). By virtue of the charge transfer mechanism, an analog voltages labelled .sup.(0) v.sub.1.sup.0. . . .sup.(0) v.sub.N.sup.b-1 are sensed at the output of each one of the b.times.N rows of the matrix array 200. To keep track of the origin of this contribution to the result, a left superscript .sup.(0) v is utilized in the notation employed herein.

In the present example, all of the products computed in the array 200 are synthesized together in accordance with corresponding powers of two to form the N elements of the transformed (output) vector. This is accomplished in an arithmetic processor illustrated in FIG. 5. The arithmetic processor is dedicated to the computation of the particular type of matrix-vector product computed in the present example. Processors other than that illustrated in FIG. 5 could be employed in making these same computations from the products produced by the array 200, and the present invention is not confined to the type of arithmetic processor illustrated in FIG. 5. Of course, in order to compute matrix-vector products different from the type computed in the present example, an arithmetic processor different from that illustrated in FIG. 5 would be employed in combination with the array 200 of FIG. 4.

At clock two cycle, the voltages sensed at each of the N rows are fed into respective pipelined A/D converters 300, each individual converter 300 having b.sub.d bits of precision (where d=log.sub.2 N, and N denotes the number of columns of A), while simultaneously A is multiplied by u.sub.1.sup.1 , . . u.sub.n.sup.1, (i.e., by the second row of charge packets in the array 230 of vector CCD/CID cells), yielding .sup.(1) v.

At clock cycle three, the digital representations of .sup.(0) v.sub.1.sup.0, are bit-mapped into a V register 305 with an appropriate bit-offset toward the most significant position. Specifically, the result or element v.sub.i.sup.l obtained during clock cycle k is offset in the appropriate V register by lk bits toward the most significant bit position (toward the leftmost bit) at clock cycle k. This offset is controlled by an offset counter 310 (which is incremented by one at each new clock cycle). Next, the voltages .sup.(1) v are fed into the A/D converters 300, and the vector u.sub.1.sup.2 , . . . u.sub.n.sup.2 multiplies A to yield .sup.(2) v.

Elements .sup.(k) v.sub.i.sup.l with same row index i are then fed into cascaded sum circuits 315 shown in FIG. 5, in parallel for all i, and pipelined over k. The cascaded sum circuits 315 are connected in a binary tree architecture of the type well-known in the art. Hence, the components v.sub.i of the product v=Au are obtained after log.sub.2 b cycles, and the overall latency is b+log.sub.2 b+3. If one needs to multiply a set of vectors u by the same matrix A, this pipelined architecture will output a new result every b clock cycles. Clearly, a far higher precision has been achieved than was previously available. An added benefit is that the refresh time overhead is significantly reduced, since in the matrix representation of FIG. 4 each electron charge packet only refers to a binary quantity.

While FIG. 4 indicates that the N.times.b.times.N cells 210 of the array 200 are formed on a single planar silicon chip, the array 200 can instead by implemented using many chips connected together, using Z-plane technology, as illustrated in FIG. 6. In such an embodiment, it is preferable to have each chip 400 assigned to a particular bit plane, in which each chip is an N-by-N CCD/CID array of CCD/CID cells 210 of the type illustrated in FIG. 1, which stores N.times.N bits or charge packets, which in the present invention, however, represent binary values only. The first chip 400-0 stores the least significant bits of all matrix elements of the N.times.N matrix, the second chip 400-1 storing the next least significant bits of all matrix elements, and so forth, and the last chip 400-(b-1) storing the most significant bits of all matrix elements. A backplane chip 405 implements the array 230 of N columns and b rows of vector CCD cells 235 of FIG. 4. The backplane chip is connected edge-wise to the column input lines of all of the bit-plane chips 400. This permits every one of the b rows of vector CCD/CID cells 235 to be input to column input lines of respective ones of the b bit-plane chips 400, greatly enhancing performance and reducing the latency of a given matrix-vector multiplication operation. The arithmetic processor of FIG. 5 could also be implemented on the backplane chip 405. The architecture of the arithmetic processor would depend upon the type of matrix-vector product to be computed. The Z-plane embodiment of FIG. 6 permits all b bits of every matrix element to be multiplied by a given vector element, and therefore is potentially much faster than the embodiment of FIGS. 4 and 5.

The pseudospectral neuralcomputing method of the present invention is now described, as it applies generally to a broad class of partial differential equations (PDE's), involving partial derivatives of various orders d of a state variable. For the sake of simplicity, the present discussion is limited to one "spatial" dimension x. Extension to multi-dimensional cases is straightforward. Thus, consider the state variable u(x,t), which is periodic over some interval which can, without loss of generality, be taken as an interval between O and L. The first step in the method is to transform the state variable u(x,t) into Fourier space with respect to x. The main advantage of this operation is that the derivatives with respect to x then become algebraic in the transformed variable. Before proceeding, the spatial interval [0,2.pi.L] is normalized to the interval [0,2.pi.]. Certain classes of PDE's can then be re-stated in terms of the new state variable v(x,t) as ##EQU1## where d is the order of the derivative of v with respect to x, v.sub.d is the partial d.sup.th derivative of v with respect to x and v.sub.t is the partial first derivative of v with respect to time. In order to numerically solve Eq. (1), the interval [0,2.pi.] is discretized by 2N equidistant points, with spacing .DELTA..sub.x =.pi./N. The function of v(x,t), which is defined only at these points, is approximated by v(x.sub.n,t), where x.sub.n =n.DELTA..sub.x, and n=0, 1 , . . . , 2N-1. The function v(x.sub.n,t) is now transformed to discrete Fourier space by ##EQU2## where k takes the values k=0, .+-.1 , . . . , .+-.N. The inversion formula is ##EQU3## This enables an efficient calculation of the derivatives of v with respect to x. In particular the d.sup.th partial derivative of v(x,t) is

where d is the order of the derivative and i=.rho.-1. In the above expression, the symbol * denotes the Schur-Hadamard product of two N.times.1 matrices. It has been shown that such a representation of derivatives of periodic functions can be related to central difference approximations of infinite order. Combining Eq. (4) with a fist order Euler approximation to the temporal derivative, yields the pseudospectral scheme for solving the PDE: ##EQU4## where .DELTA..sub.t denotes the time step size. Simple stability analysis reveals the following constraint:

In order to map the proposed numerical solution of the PDE onto a MVM architecture, Eq. (4) is expanded by substituting the appropriate terms from Eqs. (2) and (3), as follows: ##EQU5## By rearranging the terms: ##EQU6## The last term in Eq. (8) can be represented by a constant matrix which depends solely on the distance between two spatial grid points, i.e., ##EQU7## Since the function v(x.sub.n,t) is real, its derivatives should be real as well. Thus, only the real part of the matrix in Eq. (9) is needed for the computation of v.sub.m (x.sub.n,t). Hence, ##EQU8## which is readily evaluated based upon the order d of each derivative in the PDE.

The neural architecture to be employed in solving the PDE can now be specified using the foregoing results. Let each grid point represent a neuron with activation function v.sub.n (t) equal to v(x.sub.n,t). Combining Eqs. 1 and 9, the network dynamics is readily obtained: ##EQU9## Thus, the pseudo-spectral architecture consists of 2N neurons, the dynamics of which is governed by a system of coupled linear differential equations, i.e., Eq. (11). The synaptic array, T, fully interconnects all neurons, and its elements are calculated using Eq. (10).

For illustrative purposes only, a first order Euler scheme is considered for the temporal dependence. The resulting neural dynamics is then ready for implementation on the CCD/CID architecture: ##EQU10##

An architecture for implementing Equation 12 is illustrated in FIG. 7. In FIG. 7, each array 450 labelled T.sup.1, T.sup.2 , . . . etc. is preferably a CCD/CID MVM processor of the type illustrated in FIGS. 4 and 5 and stores the matrix elements defined in Equation 10. Each array 450 receives the current neuron vector v.sub.n (t) from a common register 460 labelled v.sub.n (t). Multipliers 470 multiply the outputs of the arrays 450 (T.sup.1, T.sup.2 , . . . etc.) by the appropriate coefficients from Equation 12 and an adder 480 sums the products together. A delay 485 implements the first term in Equation 12 and an adder 490 provides the final result v.sub.n (t+.DELTA..sub.t) in accordance with Equation 12. A second delay 495 provides the neuron vector v.sub.n (t) for the next computational cycle.

In order to provide a concrete framework for the proposed formalism, we focus our attention on the one-dimensional heat equation (Equation 13 below). This linear partial differential equation has the advantage of exhibiting both sufficient computational complexity, and possessing analytical solutions. This enables a rigorous benchmark of the proposed neural algorithm. The heat equation is:

where u denotes the temperature, and .alpha. is the thermal diffusivity. Partial derivatives of u with respect to time and space are denoted by u.sub.t and u.sub.x respectively. The following working example concerns an infinite slab, of thickness L, with an initial temperature distribution at time t=0 of

It is assumed that the slab is insulated at both ends, so that no heat flows through the sides. Thus, the following periodic boundary conditions apply:

Application of the pseudospectral numerical method of the invention scheme to the heat equation will now be described. This first step is to transfer u(x,t) into Fourier space with respect to x. The main advantage of this operation is that the derivatives with respect to x then become algebraic in the transformed variable. Before proceeding, the spatial period is normalized to the internal [0,2.pi.]. The scaled heat equation can then be expressed in terms of the new state variable v(x,t) as

where s is defined as 2.pi./L. In order to numerically solve Eq. (16), the interval [0,2.pi.] is discretized by 2N equidistant points, with spacing .DELTA..sub.x =.pi./N. The function v(x,t), which is defined only at these points, is approximated by v(x.sub.n,t), where x.sub.n =n.DELTA..sub.x, and n=0, 1 , . . . 2N-1. The function v(x,t) is now transformed to discrete Fourier space by: ##EQU11## where k takes the values k=0, .+-.1 , . . . .+-.N. The inversion formula is ##EQU12## This enables an efficient calculation of the derivatives of v with respect to x. In particulate

In the above expression, the symbol * denotes the Schur-Hadamard product of two N.times.1 matrices. Combining Eq. (19) with a first order Euler approximation to the temporal derivative, yields the pseudospectral scheme for solving the heat equation:

where .DELTA..sub.t denotes the time step size. The following constraints of Equation (6) applies:

In order to map the proposed numerical solution of the heat equation onto the CCD/CID architecture of FIGS. 1 or 4, the spatial derivatives are evaluated using the Fourier transform formula of Eq. (19), which is expanded by substituting the appropriate terms from Eqs. (17) and (18) to obtain ##EQU13## By rearranging the terms: ##EQU14## The last term in Eq. (23) can be represented by a constant matrix which depends solely on the distance between two spatial grid points, i.e., ##EQU15## Since the function v(x.sub.n, t) is real, its derivatives should be real as well. Thus, only the real part of the matrix in Eq. (24) is needed for the computation of v.sub.xx (x.sub.n, t). Hence, ##EQU16##

The neural architecture can now be specified for the heat equation. Let each grid point represent a neuron with activation function v.sub.n (t) equal to v(x.sub.n, t). Combining Eqs. 16 and 26, the network dynamics is readily obtained: ##EQU17## Thus, the architecture consists of 2neurons, the dynamics of which is governed by a system of coupled linear differential equations of Eq. (27). The synaptic array, T, fully interconnects all neurons, and its elements are calculated using Eq. (26).

For illustrative purposes only, a first order Euler scheme is considered for the temporal dependence. The resulting neural dynamics is then ready for implementation on the CCD/CID architecture: ##EQU18## One may further rearrange the above equation to obtain ##EQU19## denotes a synaptic matrix corresponding to the second spatial derivative in the heat equation, and includes the scaling factor, the time step, and the thermal diffusivity.

In the foregoing, the pseudo-spectral method was applied to solve linear partial differential equations such as the heat equation. The method is also applicable to quasi-linear approximations of non-linear partial differential equations such as the Korteweg-deVries equation for the soliton. The KdV or soliton equation was originally introduced to described the behavior of one-dimensional shallow water waves with small but finite amplitudes. Since its discovery, solitons have enabled many advances in areas such as plasma physics and fluid dynamics. The following description concerns the one-dimensional KdV equation, which is

where u.sub.t and u.sub.x denote partial derivatives of u with respect to time and space, respectively. If a and b are set to 6 and 1 respectively, an analytical solution to Equation 31 can be obtained for an infinite medium. This solution is:

where k and .rho..sub.0 are constants, with k>0. The above expression represents a solitary wave of amplitude 2k.sup.2 initially located at x=-.rho..sub.0 /k, moving at a velocity of 4k.sup.2. In order to numerically integrate the equation, the following boundary condition is imposed: u(x+2L,t)=u(x,t) for t in the interval [0, T] and x in the interval [-L,L].

As in the application of the pseudo-spectral method to the heat equation, the spatial interval [-L,L] is normalized to the interval [0,2.pi.] using the transformation (.pi./L)(x+L)=x.

The KdV equation can then be re-stated as:

where s=.pi./L.

Using steps analogous to those of Equations 1-11, Equation 33 becomes ##EQU20## where the matrices T.sup.1.sub.nm and T.sup.3.sub.nm are computed from Equation 10 as ##EQU21## The neural network architecture consists of 2N neurons, the dynamics of which are governed by Equation (34) defining a system of coupled non-linear differential equations. Two overlapping arrays, T.sup.1 and T.sup.3 fully interconnect all neurons.

In analogy with Equation 12, Equation 34 is recast using central difference temporal dependencies as follows: ##EQU22## and W.sup.3 is generalized for large intervals (small s) as: ##EQU23## Here, W.sup.1 and W.sup.3 are the synaptic matrices corresponding to the first and third spatial derivatives in the KdV equation, and include the scaling factor as well as the time step. Their matrix elements are individually stored in corresponding matrix cells of an MVM processor such as the CCD/CID MVM processors of FIGS. 1 or 4.

The detection of targets moving in an environment dominated by "noise" is addressed from the perspective of nonlinear dynamics. Sensor data re used to drive the Korteweg deVries (soliton) equation (Eqn. 33), inducing a resonance-type phenomenon which indicates the presence of hidden target signals. In this case, the right hand side of Eqn. 33 is not zero, but rather is set to equal the vector S.sub.n of derivatives of target sensor outputs.

Long-range detection of the motion of a target in an environment dominated by noise and clutter is a formidable challenge. Target detection problems are generally addressed from the perspective of the theory of statistical hypothesis testing. So far, existing methodologies usually fail when the signal-to-noise ratio, in dB units, becomes negative, notwithstanding the sophisticated but complex computational schemes involved. The present invention uses the phenomenology of nonlinear dynamics, not only to filter out the noise, but also to provide precise indication on the position and velocity of the target. Specifically, the invention employs the pseudo-spectral method described above to solve a driven KdV equation and achieve, by means of resonances, a dramatic enhancement of the signal to noise ratio.

To simplify the discussion, and with no loss of generality, only motion in .sup.1 space is considered. An array of N motion-detector "neurons" is fed signals, S.sub.n (t) derived from an array of N sensors. Both arrays are linear, with equally spaced elements. Let u.sub.n (t) denote the temporal response of the n-th neuron. The two overlapping synaptic arrays, W.sup.1 and W.sup.3 of Equations 39 and 40 fully interconnect the network, which obeys the following system of coupled nonlinear differential equations: ##EQU24## where W.sup.1.sub.nm and W.sup.3.sub.nm are the synaptic arrays defined in Equations 39 and 40. As in Eqn. 37, spatial organization will be considered over the interval [0.2.pi.]. The actual positions of neurons and sensors are given by the discrete values x.sub.n =n.DELTA..sub.x (n=0 , . . . N-1), with resolution .DELTA..sub.x =2.pi./N. Thus, u(x.sub.n, t) is written as u.sub.n (t). For convenience N is even. The sensor derivative inputs to the network are denoted by S.sub.n (t).

The corresponding hardware architecture is illustrated in FIG. 8. Synaptic arrays 500, 510 are each integrated circuits (chips) embodying CCD/CID arrays of the type illustrated in FIG. 4 and store the matrix elements of W.sup.1 and W.sup.3 respectively. Arithmetic processors 515, 520 are each chips embodying arithmetic processors of the type illustrated in FIG. 5. Registers 525 and 530 hold the current values of u(t). A detector array 535 furnishes the derivative signal vector S.sub.n (t). The second term of Equation 41 is implemented by the CCD/CID array 500, the arithmetic processor 515 and a multiplier 540. The third term of Equation 41 is implemented by the CCD/CID array 510, the arithmetic processor 520 and an adder 545. The last term of Equation 41 is implemented by the detector array 535 and an adder 550. A delay 555 implements the first term of Equation 41.

The results of a computer simulation of FIG. 8 are illustrated in FIGS. 9A, 9B, 9C, 9D, 10A, 10B, 10C, 10D, 11A, 11B, 11C and 11D. Each drawing of FIGS. 9A, 9B, 9C, 9D, 10A, 10B, 10C, 10D, 11A, 11B, 11C and 11D includes four plots. The lower left corner (FIGS. 9A, 10A, 11A) shows the sensors input to the network. Since data are simulated, the contributions from the target and background noise are plotted separately, for illustrative purposes, even though the network actually receives only their combined value. The Fourier spectrum of the total signal is given in the upper left corner (FIGS. 9B, 10B and 11B). The network's output is plotted in the lower right corner (FIGS. 9C, 10C, 11C). Finally, the solution in absence of noise is shown in the upper right corner (FIGS. 9D, 10D, 11D).

The results are plotted after 10 time steps (.DELTA..sub.t =0.005), and were obtained using a 64-element sensor array. FIG. 9 indicates that when only noise is fed to the network, no spurious results emerges. FIG. 10 illustrates the spectral network's detection capability for a target moving in a noise and clutter background characterized by a per pulse signal-to-noise (SNR) ratio of approximately 0 dB. Conventional detection methods usually reach their breaking point in the neighborhood of such SNR rations. Finally, FIG. 11 presents results for a case where the SNR drops below -10 dB. The above SNR is by no way the limit, as can be inferred from the still excellent quality of the detection peak. Furthermore, multiple layers of the spectral neural architecture provide a "space-time" tradeoff for additional enhancement. The following presents a possible physical explanation of the detection phenomenon.

In interpreting the observed detection phenomenon of Eqn. 41, it is assumed that the sensor input to the network consists of two parts: (1) .rho.(x,t), space-dependent random oscillations; and (2) a "target" .THETA.(x,t), the value of which is nonzero only over a few "pixels" of the sensor array. Since the KdV equations are non dissipative, and energy is pumped into the system via the S(x,t) term, active stems have to be taken to avoid unbounded growth of the solutions. Furthermore, if v denotes the target's velocity, then

Thus, based upon the properties of the homogeneous KdV equation (Eqn. 34), the target signal .THETA.(x,t) "resonates" with the "eigen-solitons" of the homogeneous equation. Hence, it will be amplified, while the random components .rho.(x,t) will be dispersed. Furthermore, such a "resonance" should not depend on the target velocity, since the velocity of the "eigen-solitons" of the KdV equations are not prescribed. In other words, the proposed methodology can detect targets over a wide range of velocities.

While the invention has been described in detail by specific reference to preferred embodiments, it is understood that variations and modifications thereof may be made without departing from the true spirit and scope of the invention.