Signal processing of multi-channel data

[0001] The present invention relates to signal processing, and is more particularly related to linear prediction.

[0002] Signals can represent information from any source that generates data, relating to electromagnetic energy to stock prices. Analysis of these signals is the focus of signal processing theory and practice. Linear prediction is an important signal processing technique that provides a number of capabilities: (1) prediction of the future of a signal from its past; (2) extraction of important features of a signal; and (3) compression of signals. The economic value of linear prediction is incalculable as its prevalence in industry is enormous.

[0003] It is observed that many important signals are “multi-channel” in that the signals are gathered from many independent sources; e.g., time series. For example, multi-channel data stem from the process of searching for oil, which requires measuring the earth at many locations simultaneously. Also, measuring the motions of walking (i.e., gait) requires simultaneously capturing the positions of many joints. Further, in a video system, a video signal is a recording of the color of every pixel on the screen at the same moment; essentially each pixel is essentially a separate “channel” of information. Linear prediction can be applied to all of the above disparate applications.

[0004] Conventional linear prediction techniques have been inadequate in the treatment of multi-channel time series, particularly, when the dimensionality is in the order is above three. There are traditional approaches of linear prediction for multi-channel signals, but are not effective in addressing the technical difficulties that are caused by the interactions of the sources of data. In single source signals, such as like voice, these difficulties are not encountered. The conventional techniques assume that the autocorrelation matrix of the data is invertible or can be made invertible by simple methods, which is rarely valid for real multi-channel data.

[0005] Also, such traditional approaches do not use the structural information available through modeling multi-dimensional geometry in a more sophisticated manner than merely as arrays of numbers. In addition, these approaches fail to take into account the phenomenon of time warping, which, for example, is critical to successful modeling of biometric time series. Further, conventional linear prediction techniques are based on a statistical foundation for linear prediction, which is not well suited for motion, video and other types of multi-channel data.

[0006] Further, it is recognized that most real multi-channel data are highly correlated. Under the conventional approaches, the popular linear prediction algorithm, known as the Levinson algorithm, cannot be applied to highly correlated channels.

[0007] Therefore, there is a need to provide a framework for extending applicability of linear prediction techniques. Additionally, there is a need for an approach to predict/compress/encrypt multi-channel multi-dimensional time series, particularly series with high correlation.

[0008] These and other needs are addressed by the present invention in which non-commutative approaches to signal processing are provided. In one embodiment, quaternions are used to represent multi-dimensional data (e.g., three- and four-dimensional data, etc.). Additionally, an embodiment of the present invention provides a linear predictive coding scheme (e.g., based on the Levinson algorithm) that can be applied to a wide class of signals in which the autocorrelation matrices are not invertible and in which the underlying arithmetic is not commutative. That is, the linear predictive coding scheme can handle singular autocorrelations, both in the commutative and non-commutative cases. Random path modules are utilized to replace the statistical basis of linear prediction. The present invention, according to one embodiment, advantageously provides an effective approach for linearly predicting multi-channel data that is highly correlated. The approach also has the advantage of solving the problem of time-warping.

[0009] In one aspect of the present invention, a method for providing linear prediction is disclosed. The method includes collecting multi-channel data from a plurality of independent sources, and representing the multi-channel data as vectors of quaternions. The method also includes generating an autocorrelation matrix corresponding to the quaternions. The method further includes outputting linear prediction coefficients based upon the autocorrelation matrix, wherein the linear prediction coefficients represent a compression of the collected multi-channel data.

[0010] In another aspect of the present invention, a method for supporting video compression is disclosed. The method includes collecting time series video signals as multi-channel data, wherein the multi-channel data is represented as vectors of quaternions. The method also includes generating an autocorrelation matrix corresponding to the quaternions, and outputting linear prediction coefficients based upon the autocorrelation matrix.

[0011] In another aspect of the present invention, a method of signal processing is provided. The method includes receiving multi-channel data, representing multi-channel data as vectors of quaternions, and performing linear prediction based on the quaternions.

[0012] In another aspect of the present invention, a method of performing linear prediction is provided. The method includes representing multi-channel data as a pseudo-invertible matrix, generating a pseudo-inverse of the matrix, and outputting a plurality of linear prediction weight values and associated residual values based on the generating step.

[0013] In another aspect of the present invention, a computer-readable medium carrying one or more sequences of one or more instructions for performing signal processing is disclosed. The one or more sequences of one or more instructions include instructions which, when executed by one or more processors, cause the one or more processors to perform the steps of receiving multi-channel data, representing multi-channel data as vectors of quaternions, and performing linear prediction based on the quaternions.

[0014] In yet another aspect of the present invention, a computer-readable medium carrying one or more sequences of one or more instructions for performing signal processing is disclosed. The one or more sequences of one or more instructions include instructions which, when executed by one or more processors, cause the one or more processors to perform the steps of representing multi-channel data as a pseudo-invertible matrix, generating a pseudo-inverse of the matrix, and outputting a plurality of linear prediction weight values and associated residual values based on the generating step.

[0015] Still other aspects, features, and advantages of the present invention are readily apparent from the following detailed description, simply by illustrating a number of particular embodiments and implementations, including the best mode contemplated for carrying out the present invention. The present invention is also capable of other and different embodiments, and its several details can be modified in various obvious respects, all without departing from the spirit and scope of the present invention. Accordingly, the drawing and description are to be regarded as illustrative in nature, and not as restrictive.

[0022] A system, method, and software for processing multi-channel data by non-commutative linear prediction are described. In the following description, for the purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the present invention. It is apparent, however, to one skilled in the art that the present invention may be practiced without these specific details or with an equivalent arrangement. In other instances, well-known structures and devices are shown in block diagram form in order to avoid unnecessarily obscuring the present invention.

[0023] The present invention has applicability to a wide range of fields in which multi-channel data exist, including, for example, virtual reality, doppler radar, voice analysis, geophysics, mechanical vibration analysis, materials science, robotics, locomotion, biometrics, surveillance, detection, discrimination, tracking, video, optical design, and heart modeling.

[0024] FIG. 1 is a diagram of a system for providing linear prediction, according to an embodiment of the present invention. As shown in FIG. 1, a multi-channel data source 101 provides data that is converted to quaternions by a data representation module 103. Quaternions have not been employed in signal processing, as conventional linear prediction techniques cannot process quaternions in that these techniques employ the concept of numbers, not points. According to one embodiment of the present invention, quaternions can be parsed into a rotational part and a scaling part; this construct, for example, can correct time warping, as will be more fully described below.

[0025] These quaternions are then supplied to a non-commutative linear predictor 105, which generates the linear prediction matrix 107 of weights and associated residuals. The linear predictor 105, in an exemplary embodiment, provides a generalization of the Levinson algorithm to process non-invertible autocorrelation matrices over any ring that admits compact projections. Linear predictive techniques conventionally have been presented in a statistical context, which excludes the majority of multi-channel data sources to which the linear predictor 105 is targeted.

[0026] The signal processing of spatial time series has been traditionally limited by the lack of a sophisticated link between the signal processing algebra and the spatial geometry. The ordinary algebra of the real or complex numbers satisfies the commutative law a×b=b×a and the law of inverses: for every non-zero number a there is a number

[0027] $\frac{1}{a}$
for which

[0028] $a\times \left(\frac{1}{a}\right)=\left(\frac{1}{a}\right)\times a=1.$
However, these properties fail for the quaternions and for three-dimensional multi-channel signal processing. The theories of hermitian regular rings and compact projections allow important signal processing techniques to be utilized in such situations.

[0029] One of the major application areas of the invention is to video image processing. To enable this application, color data needs to be correctly represented as four-dimensional spatial points. Photopic coordinates are four-dimensional analogs of the common RGB (Red-Green-Blue) colormetric coordinates.

[0030] Also, in gait analysis, for example, each joint reports where it currently is located. In the oil exploration example, each of many sensors spread over the area that is being searched sends back information about where the surface on which it is sitting is located after the geologist has set off a nearby explosion. The cardiology example requires knowing, for many structures inside and around the heart, how these structures move as the heart beats.

[0031] Even the video example can be seen that way because each pixel on the screen is reporting its color at every moment of time. However, a “color” is not a simple number: it is actually (at least) 3 numbers such as the amount of red, blue, and green (RGB) light needed to make that color. Those three numbers are usually thought of as being in a “color space” which is a kind of abstract space like three-dimensional space.

[0032] As mentioned, the present invention, according to one embodiment, represents each such point in space by a mathematical object called a “quaternion.” Quaternions can describe special information, such as rotations, perspective drawing, and other simple concepts of geometry. If a signal, such as the position of a joint during a walk is described using quaternions, it reveals structure in the signal that is hidden such as how the rotation of the knee is related to the rotation of the ankle as the walk proceeds.

[0033] FIGS. 2A and 2B are diagrams of multi-channel data capable of being processed by the system of FIG. 1. As shown in FIG. 2A, many practical datasets comprise time series . . . x_n−2, x_n−1, x_n of data vectors where, at each time n, the datum x_n is a vector

[0034] $x_n=\left(\begin{array}{c}{x}_n\left(1\right)\\ x_n\left(2\right)\\ ⋮\\ x_n\left(K\right)\end{array}\right)$
of three-dimensional measurements. Each component x_n(k) represents the measurement of a single channel and is itself composed of three separate real numbers x_n(k)=(x_n(k)¹x_n(k)²x_n(k)³) corresponding to the three dimensions of whatever system that is being measured.

[0035] It is clear that cross-channel measurements can be represented as a list, x_n:

[0036] $x_n=\left(\begin{array}{ccc}\left(\begin{array}{c}{x_n\left(1\right)}^1\\ {x_n\left(2\right)}^1\\ ⋮\\ {x_n\left(K\right)}^1\end{array}\right)& \left(\begin{array}{c}{x_n\left(1\right)}^2\\ {x_n\left(2\right)}^2\\ ⋮\\ {x_n\left(K\right)}^2\end{array}\right)& \left(\begin{array}{c}{x_n\left(1\right)}^3\\ {x_n\left(2\right)}^3\\ ⋮\\ {x_n\left(K\right)}^3\end{array}\right)\end{array}\right),$
such as the RGB bitplanes of video and, in fact, this is usually how three-dimensional datasets are generated. However, the former representation is conceptually more basic.

[0037] As seen in FIG. 2B, a time series relating to the prices of stocks, for example, exist, and can be viewed as a single multi-channel data. In this example, three sources 201, 203, 205 can be constructed as a single vector based on time, t.

[0038] According to one embodiment of the present invention, multi-channel can be represented as quaternions. Specifically, the present invention provides an approach for analyzing and coding such time series by representing each measurement x_n(j) using the mathematical construction called a quaternion.

[0039] FIG. 3 is a flow chart of a process for representing multi-channel data as quaternions, according to an embodiment of the present invention. In step 301, multi-channel data is collected and then represented as quaternions, as in step 303. These quaternions, per step 305, are then output to a linear predictor (e.g., predictor 105 of FIG. 1).

[0040] As used herein, the quaternion algebra is denoted . Quaternions are four-dimensional generalizations of the complex numbers and may be viewed as a pair of complex numbers (as well as many other representations). Quaternions also have the standard three-dimensional dot- and cross-products built into their algebraic structure along with four-dimensional vector addition, scalar multiplication, and complex arithmetic.

[0041] The quaternions have the arithmetical operations of +,−,×, and ÷ for non-0 denominators defined on them and so provide a scalar structure over which vectors, matrices, and the like may be constructed. However, the peculiarity of quaternions is that multiplication is not commutative: in general, q×r≠r×q for quaternions q,r and thus forms a division ring, not a field.

[0042] The present invention, according to one embodiment, presented herein stems from the observation that many traditional signal processing algorithms, especially those pertaining to linear prediction and linear predictive coding, do not depend on the commutative law holding among the scalars once these algorithms are carefully analyzed to keep track of which side (left or right) scalar multiplication takes place.

[0043] As a result, a three- (or four-) dimensional data point can be thought of as a single arithmetical entity rather than a list of numbers. There are great advantages to be gained, both conceptually and practically, by doing so.

[0044] As mentioned previously, the application of present invention spans a number of disciplines, from biometrics to virtual reality. For instance, all human control devices from the mouse or gaming joystick up to the most complex virtual reality “suit” are mechanisms for translating spatial motion into numerical time series. One example is a “virtual reality” glove that contains 22 angle-sensitive sensors arrayed on a glove. Position records are sent from the glove to a server at 150 records/sensor/sec at the RS-232 rate of 115.2 kbaud. After conversion to rectangular coordinates, this is precisely a 22-channel time series . . . x_n−2, x_n−1, x_n,

[0045] $x_n=\left(\begin{array}{c}{x}_n\left(1\right)\\ x_n\left(2\right)\\ ⋮\\ x_n\left(22\right)\end{array}\right)$
of three-dimensional data as discussed above.

[0046] The high data rate and sensor sensitivity of the virtual glove is sufficient to characterize hand positions and velocities for ordinary motion. However, the human hand is capable of “extraordinary” motion; e.g., a skilled musician or artisan at work. For example, both pianists and painters have the concept of “touch”, an indefinable relation of the hand/finger system to the working material and which, to the trained ear or eye, characterizes the artist as well as a photograph or fingerprint. It is just such subtle motions, which unerringly distinguish human actions from robotic actions.

[0047] Even to begin the modeling and reproduction of the true human hand, much higher data rates, much more precise sensors, and much denser sensor array are required. The numbers are comparable, in fact, to the data rates, volume, and density of the nervous system connecting the hand to the brain. At such levels, efficient storing and transmission of such multi-channel data become critical. It is not sufficient to save bandwidth by transmitting only every tenth or hundredth hand position of a pilot landing a jet fighter on the flight deck of a carrier. Instead, the time series need to be globally compressed so that actual redundancy (introduced by inertia and physiological/geometric constraints) but not critical information is removed.

[0048] Multi-channel analysis is also utilized in geophysics. Geophysical explorers, like special effects people in cinema, are in the enviable position of being able to set off large explosions in the course of their daily work. This is a basic mode of gathering geophysical data, which arrives from these earth-shaking events (naturally occurring or otherwise) in the form of multi-channel time series recording the response of the earth's surface to the explosions. Each channel represents the measurements of one sensor out of a strategically-designed array of sensors spread over a target area.

[0049] While the input data series of any one channel is typically one-dimensional, representing the normal surface strain at a point, the target series is three-dimensional; namely, the displacement vector of each point in a volume. Geophysics is, more than most sciences, concerned with inverse problems: given the boundary response of a mechanical system to a stimulus, determine the response of the three-dimensional internal structure. As oil and other naturally occurring resources become harder to find, it is imperative to improve the three-dimensional signal processing techniques available.

[0050] Similar to geophysicists, mechanical engineers examine system response measurements. Typically, a body is covered in a multi-channel network of strain or motion sensors and shakers is attached at selected points. The data usually is transferred to a finite-element model of the system, which is a triangularization of the three-dimensional physical system. Abstractly, these finite-element datasets are nothing more than the multi-channel three-dimensional time series.

[0051] Multi-channel analysis also has applicability to biophysics. If a grid is placed over selected points of photographed animals' bodies, and concentrated especially around the joints, time series of multi-channel three-dimensional measurements can be generated from these historical datasets by standard photogrammetric techniques.

[0052] The human knee is a complex mechanical system with many degrees of freedom most of which are exercised during even a simple stroll. This applies to an even greater degree to the human spine, with its elegant S-shape, perfectly designed to carry not only the unnatural upright stance of homo sapiens but to act as a complex linear/torsional spring with infinitely many modes of behavior as the body walks, jumps, runs, sleeps, climbs, and, not least of all, reproduces itself. Many well-known neurological diseases, such as multiple sclerosis, can be diagnosed by the trained diagnostician simply by visual observation of the patient's gait.

[0053] Paleoanthropologists use computer reconstructions of hominid gaits as a basic tool of their trade, both as an end product of research and a means of dating skeletons by the modernity of the walk they support. Animators are preeminent gait modelers, especially these days when true-to-life non-existent creatures have become the norm.

[0054] The present invention also applicability to biometric identification. Closely related to the previous example is the analysis of real human individuals' walking characteristics. It is observed that people frequently can be identified quite easily at considerable distances simply by their gait, which seems as characteristic of a person as his fingerprints. This creates some remarkable possibilities for the identification and surveillance of individuals by extracting gait parameters as a signature.

[0055] It might be possible, for example, to establish the identity of a criminal suspect through analysis of gait characteristics from closed circuit television (CCTV) recording, even when the quality of these videos is too poor to isolate facial structure. A system could be constructed that would follow a particular individual through, say, a crowded airport or cityscape by identifying his walking signature via CCTV. An ordinary disguise, of course, will not fool such a system. Even the conscious attempt to walk differently may not succeed because the primary determinants of gait (such as the particular mechanical properties of the spine/pelvis interface) may be beyond conscious control.

[0056] The present invention, additionally, is applicable to detection, discrimination, and tracking of targets. There are many targets which move in three spatial dimensions and which it may be desirable to detect and track. For example, a particular aircraft or an enemy submarine in the ocean. Although there are far fewer channels than in gait analysis, these target tracking problems have a much higher noise floor.

[0057] There are many well-known techniques of adapting linear prediction to noisy signals, one of the simplest yet most effective being to manually adjust the diagonal coefficients of the autocorrelation matrix.

[0058] Multi-channel analysis can also be applied to video processing. Spatial measurements are not the only three-dimensional data which has to be compressed, processed, and transmitted. Color is (in the usual formulations) inherently three-dimensional in that a color is determined by three values: RGB, YUV (Luminance-Bandwidth-Chrominance), or any of the other color-space systems in use.

[0059] A video stream can be modeled by the same time series . . . x_n−2, x_n−1, x_n approach that has been traditionally employed, except that now a channel corresponds to a single pixel on the viewing screen:

[0060] $x_n=\left(\begin{array}{ccc}{C}_n\left(11\right)& \cdots & C_n\left(1N\right)\\ ⋮& ⋰& ⋮\\ C_n\left(\mathrm{M1}\right)& \cdots & C_n\left(\mathrm{MN}\right)\end{array}\right)$

[0061] where C_n(jk)=(C_n(jk)^R C_n(jk)^G C_n(jk)^B) are the three color coordinates at time n in, for example, the RGB system of pixel j,k out of a total resolution of (M×N) pixels.

[0062] As mentioned previously, many hardware systems require the data to be arranged in the dual form of three value planes rather than planes of three values. With the large quantity of data represented by . . . x_n−2, x_n−1, x_n, compression is the key to successful video manipulation. For example, there is increasing pressure for corporate intranets to carry internal video signals and, for these applications, security is a critical necessity from the outset.

[0063] According to one embodiment, the present invention introduces the concept of photopic coordinates; it is shown that, just as in spatial data, color data is modeled effectively by quaternions. This construct permits application of the non-commutative methods to color images and video a reanalysis of the usual color space has to be performed, recognizing that color space has an inherent four-dimensional quality, in spite of the three-dimensional RGB and similar systems.

[0064] Many signal processing problems are presented in the form of overlapping frames laid over a basic single-channel time series:

[0065] $\begin{array}{cccc}{x}_1& x_2& \cdots & x_K\end{array}{x}_{K+1}\cdots x_n\cdots$ $x_1\cdots \begin{array}{cccc}{x}_{d+1}& x_{d+2}& \cdots & x_{d+K}\end{array}\cdots x_n\cdots$ $x_1{x}_2⋮.x_1{x}_2\cdots \begin{array}{cccc}{x}_{\mathrm{md}+1}& x_{\mathrm{md}+2}& \cdots & x_{\mathrm{md}+K}\end{array}\cdots$ $⋮$

[0066] High-resolution spectral analysis by linear prediction or some other method is performed separately within each frame

[0067] $\begin{array}{cccc}{x}_{\mathrm{md}+1}& x_{\mathrm{md}+2}& \cdots & x_{\mathrm{md}+K}\end{array}$
and then the resulting power spectra P₀(ω), P₁(ω), . . . , P_m(ω), . . . are analyzed as a new data sequence.

[0068] This is the traditional approach in voice analysis where the resulting spectra are presented in the well-known spectrogram form. However, it is used in many other applications such as the Doppler radar analysis of rotating bodies in which the distances of reflectors from the axis of rotation can be deduced from the instantaneous spectra of the returned signal.

[0069] More generally, this frame-based spectral analysis can be regarded as the demodulation of an FM (Frequency Modulation) signal because the information that is to be extracted is contained in the instantaneous spectra of the signal. Unfortunately, this within-frame approach ignores some of the most important information available; namely the between-frame correlations.

[0070] For example, in the rotating Doppler radar problem, a single rotating reflector gives rise to a sinusoidally oscillating frequency spike in the spectra sequence P₀(ω), P₁(ω), . . . , P_m(ω), . . . . The period of oscillation of this spike is the period of rotation of the reflector in space while the amplitude of the spike's oscillation is directly proportional to the distance of the reflector from the axis of rotation. These oscillation parameters cannot be read directly from any individual spectrum P_m(ω) because they are properties of the mutual correlations between the entire sequence P₀(ω), P₁(ω), . . . , P_m(ω), . . . .

[0071] This point is brought out especially well in the presence of noise which, as is well-known, has a strongly deleterious effect on any high-resolution spectral analysis method. An individual spectrum P_m(ω) may not exhibit any discernable spike but since it is known that there is an underlying oscillation in the series P₀(ω), P₁(ω), . . . , P_m(ω), . . . , a way exists to combine these spectra to filter out the cross-frame noise.

[0072] It is recognized that by imposing the frame structure on the time sequence, the signal is transformed into a multi-channel sequence:

[0073] $\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_K\end{array},\begin{array}{c}{x}_{d+1}\\ x_{d+2}\\ ⋮\\ x_{d+K}\end{array},\cdots ,\begin{array}{c}{x}_{\mathrm{md}+1}\\ x_{\mathrm{md}+2}\\ ⋮\\ x_{\mathrm{md}+K}\end{array},\cdots$
with the number of channels K equal to the frame width.

[0074] As is more fully described below, linear predictive analysis of such a multi-channel sequence gives rise to coefficients a₁, . . . , a_m, . . . which are (K×K) matrices rather than single scalars. Thus, the spectra P_m(ω) produced by these coefficients are themselves (K×K) matrices.

[0075] However, the correlations that are sought after, such as the oscillation patterns produced by rotating radar reflectors, cause these power spectra matrix sequences P₀(ω), P₁(ω), . . . , P_m(ω), . . . to become singular; i.e., the autocorrelation matrices of P₀(ω), P₁(ω), . . . , P_m(ω), . . . (which are matrices whose entries are themselves matrices) becomes non-invertible. In fact, the non-inevitability of this matrix is equivalent to cross-spectral correlation.

[0076] Unfortunately, the prior approaches to linear prediction break down at this exact point because these conventional approaches cannot handle the problem of channel degeneracy.

[0077] The present invention, according to one embodiment, advantageously operates in the presence of highly degenerate data.

[0078] As noted, the present invention can be utilized in the area of optics. It has been understood that optical processing is a form of linear filtering in which the two-dimensional spatial Fourier transforms of the input images are altered by wavenumber-dependent amplitudes of the lens and other transmission media. At the same time, light itself has a temporal frequency parameter ν which determines the propagation speed and the direction of the wave fronts by means of the frequency-dependent refractive index. Thus, the abstract optical design and analysis problem is determining the relation between the four-component wavevector ({right arrow over (σ)},ν) and the on the four-component space-time vector ({right arrow over (x)},t) on each point of a wavefront as it moves through the optical system.

[0079] Both ({right arrow over (σ)},ν) and ({right arrow over (x)},t) for a single point on a wavefront can be viewed as series of four-dimensional data, and thus, a mesh of points on a wavefront generates two sets of two-dimensional arrays of four-dimensional data. As is seen, ({right arrow over (σ)},ν),({right arrow over (x)},t) are naturally structured as quaternions. There are many possibilities for joint linear predictive analysis of these series. In particular, estimating the four-dimensional power spectra by solving for the all-pole filter produced by the linear prediction model.

[0080] Passing from two-dimensional arrays of three-dimensional data, there are many applications which require three-dimensional arrays of three-dimension data. For example, the stress of a body is characterized by giving, for every point (x,y,z) inside the unstressed material, the point (x+εx,y+εy,z+εy) to which (x,y,z) has been moved. If a uniform grid of points (lΔx,mΔy,nΔz), {l,m,n}⊂³ defines the body, then the three-dimensional array

[0081] $\left(\begin{array}{ccc}& ⋮& \\ \cdots & {\left(\delta x,\delta y,\delta z\right)}_{l,m,n}& \cdots \\ & ⋮& \end{array}\right)\to )$

[0082] of three-dimensional data approximates the stress. For example, from this matrix, an approximation of the stress tensor may be derived.

[0083] A good example of the use of these ideas is three-dimensional, dynamic modeling of the heart. The stress matrix can be obtained from real-time tomography and then linear predictive modeling can be applied. This has many interesting diagnostic applications, comparable to a kind of spatial EKG (Electrocardiogram).

[0084] As is discussed later, the system response of the quaternion linear filter is a function of two complex values (rather than one as in the commutative situation). Thus the “poles” of the system response really is a collection of polar surfaces in ×≅⁴. Because of the strong quasi-periodicities in heart motion and because the linear prediction filter is all-pole, these polar surfaces can be near to the unit 3-sphere (the four-dimensional version of the unit circle) in ⁴.

[0085] The stability of the filter is determined by the geometry of these surfaces, especially by how close they approach the 3-sphere. It is likely that this can be translated into information about the stability of the heart motion, which is of great interest to cardiologists.

[0086] FIG. 4 is a flowchart of the operation for performing non-commutative linear prediction in the system of FIG. 1. Linear prediction (LP) has been a mainstay of signal processing, and provides, among other advantages, compression and encryption of data. Linear prediction and linear predictive coding, according to one embodiment of the present invention, requires computation of an autocorrelation matrix of the multi-channel data, as in step 301. While theoretically creating the possibility of significant compression of multi-channel sets, such high degrees of correlation also create algorithmic problems because it causes the key matrices inside the algorithms to become singular or, at least, highly unstable. This phenomenon can be termed “degeneracy” because it is the same effect which occurs in many physical situations in which energy levels coalesce due to loss of dimensionality.

[0087] Degeneracy cannot be removed simply by looking for “bad” channels and eliminating them. For one thing, such a scheme is too costly in time, and fundamentally flawed, because degeneracy is a global or system-wide phenomenon. The problem of degeneracy of multi-channel data has generally been ignored by algorithm designers. For example, traditional approaches only consider the case in which the autocorrelation matrices are either non-singular (another way of saying the system is not degenerate) or that the singularity can be confined to a few deterministic channels. Without this assumption, the popular linear prediction method, referred to as the Levinson algorithm, fails in its usual formulation.

[0088] Real multi-channel data, as discussed above, can be expected to be highly degenerate. The present invention, according to one embodiment, can be used to formulate a version of the Levinson algorithm that does not assume non-degenerate data. This is accomplished by examining the manner in which matrix inverses enter into the algorithm; such inverses can be replaced by pseudo-inverses. This is an important advance in multi-channel linear prediction even in the standard commutative scalar formulations.

[0089] In step 303, pseudo-inverses of the autocorrelation matrix are generated, thereby overcoming any limitations stemming for the non-inevitability problem. The linear predictor then outputs the linear prediction matrix containing the LP coefficients and residuals, per step 305.

[0090] The general idea of compression is that any data set contains hidden redundancy which can be removed, thus reducing the bandwidth required for the data's storage and transmission. In particular, predictive coding removes the redundancy of a time series . . . x_n−2, x_n−1, x_n by determining a predictor function ( ) and a new residual data series . . . e_n−2, e_n−1, e_n for which
x_n=(x_n−1,x_n−2, . . . )+e_n
for every n in an appropriate range. Ideally, ( ) will depend on relatively few parameters, analogous to the coefficients of a system of differential equations and which are transmitted at the full bit-width, while . . . e_n−2, e_n−1, e_n will have relatively low dynamic range and thus can be transmitted with fewer bits/symbol/time than the original series. The series, . . . e_n−2, e_n−1, e_n, can be thought of as equivalent to the series . . . x_n−2, x_n−1, x_n but with the deterministic redundancy removed by the predictor function ( ). Equivalently, . . . e_n−2, e_n−1, e_n is “whiter” than . . . x_n−2, x_n−1, x_n; i.e., has higher entropy per symbol.

[0091] The compression can be increased by allowing lossy reconstruction in which only a fraction (possibly none) of the residual series . . . e_n−2, e_n−1, e_n is transmitted/stored. The missing residuals are reconstructed as 0 or some other appropriate value.

[0092] Encryption is closely associated with compression. Encryption can be combined with compression by encrypting the ( ) parameters, the residuals . . . e_n−2, e_n−1, e_n, or both. This can be viewed as adding encoded redundancy back into the compressed signal, analogous to the way error-checking adds unencoded redundancy.

[0093] Linear prediction and linear predictive coding use a finite linear function
(x_n−1,x_n−2,x_n−3, . . . )=−a₁x_n−1−a₂x_n−2−a₃x_n−3 . . . −a_Mx_n−M
with constant coefficients as the predictor.

[0094] So defining a₀=1, the full LP model of order M is

[0095] $\sum _{m=0}^M{a}_m{x}_{n-m}=e_n$

[0096] It is noted that when each x_n is a K-channel datum, the coefficients a_m must be (K×K) matrices over the scalars (typically ,, or ).

[0097] A number of non-LP coding schemes exists, such as the Fourier-based JPEG (Joint Photographic Experts Group) standard. The LP models have a universality and tractability which make them benchmarks.

[0098] Linear prediction becomes statistical when a probabilistic model is assumed for the residual series, the most common being independence between times and multi-normal within a time; that is, between channels at a single moment of time when each x_n is a multi-channel data sample.

[0099] The property enjoyed by the multi-normal density

[0100] $\varphi \left(x_1,\cdots ,x_n\right)=\varphi \left(\overrightarrow{x}\right)=\frac{1}{{\left(2\pi \right)}^{n/2}}\frac{1}{\sqrt{\det \Sigma }}{e}^{-\frac{1}{2}{\left(\overrightarrow{x}-\overrightarrow{\mu }\right)}^T{\Sigma }^{-1}\left(\overrightarrow{x}-\overrightarrow{\mu }\right)},$
where Σ is the covariance matrix and {right arrow over (μ)} the mean of {right arrow over (x)}, and no other distribution is that uncorrelated multi-normal random variables are statistically independent. As a result, “independent” in the sense of linear algebra is identical to “independent” in the sense of probability theory. By linearly transforming the variables to the principal axes determined by the eigenstructure of Σ, consideration can be narrowed to independent, normally distributed random variables. The residuals can be tested for significance using standard χ²- or F-tests, analysis of variance (ANOVA) tables can be constructed, and the rest.

[0101] In essence, then, any advancement of linear predictive coding must either improve the linear algebra or improve the statistics or both.

[0102] The present invention advances the linear algebra by introducing non-commutative methods, with the quaternion ring as a special case, into the science of data coding. The present invention also advances the statistics by reanalyzing the basic assumptions relating linear models to stationary, ergodic processes. In particular, it is demonstrated by analyzing source texts that linear prediction is not a fundamentally statistical technique and is, rather, a method for extracting structured information from structured messages.

[0103] Like all signal processing methodologies, the three-dimensional, non-commutative technique is a series of modeling “choices,” not just one algorithm applicable to all situations. As a result of this and due to the unfamiliarity of many of the mathematical concepts being used, an attempt is made to provide a reasonably self-contained presentation of the context in which the modeling takes place.

[0104] In statistical signal processing, LP appears as autoregressive models (AR). These are a special case of autoregressive-moving average models (ARMA) which, unlike AR models, have both poles and zeros; i.e. modes and anti-modes. For example, in radar applications, the same general class of techniques are usually called autoregressive spectral analysis and have found diverse applications including target identification through LP analysis of Doppler shifts.

[0105] As pointed out previously, the K-channel linear predictive model is as follows:

[0106] $\underset{m=0}{\sum ^M}{a}_m{x}_{n-m}=e_n,$
which requires the coefficients a_m to be (K×K) matrices which, in general, do not commute: a·b=b·a. As is discussed below, when the entries of the matrices a_m themselves are commutative, the non-commutativity of the a_m can be controlled at the determinants since det(a·b)=det(b·a) even when a·b=b·a.

[0107] However, once the matrices are composed of non-commutative entries, the determinant is no longer useful. This results, for example, if higher-order prediction is to be performed in which multiple channels of series (which are themselves multi-channel series are utilized). This is not an abstraction: many real series are presented in this form. For example, it may be the case that the multi-channel readings of geophysical experiments from many separate locations are used and it is desired to assemble them all into a single predictive model for, say, plate tectonic research. It is not the case that the model derived by representing all channels into a large, flat matrix is the same as that obtained by regarding the coefficients a_m as matrices whose entries are also matrices.

[0108] The general linear prediction problem is thus concerned with the algebraic properties of the set (n,m,A) of (n×m) matrices whose entries are in some scalar structure A. Appropriate scalar structures is discussed in below with respect to quaternion representations. In many cases, however, A is itself a matrix structure (k,l,B). There is thus a tendency to regard aε(n,m,A), with A=(k,l,B), as “really” structured as aε(nk,ml,B):

[0109] $\begin{array}{c}\underset{\downarrow }{\overset{\uparrow }{n}}\stackrel{\leftarrow m\to }{\left(\begin{array}{ccc}{a}_{11}& \cdots & a_{1m}\\ ⋮& ⋰& ⋮\\ a_{n1}& \cdots & a_{nm}\end{array}\right)},a_{v\mu }=\underset{\downarrow }{\overset{\uparrow }{k}}\stackrel{\leftarrow l\to }{\left(\begin{array}{ccc}{a}_{v\mu ,11}& & a_{v\mu ,1l}\\ & & \\ & & \\ a_{v\mu ,k1}& & a_{v\mu ,\mathrm{kl}}\end{array}\right)}\\ \Updownarrow \\ \underset{\downarrow }{\overset{\uparrow }{\mathrm{nk}}}\stackrel{\leftarrow ml\to }{\left(\begin{array}{cccccccccc}{a}_{11,11}& \cdots & a_{12,11}& \cdots & & & & & \cdots & a_{1m,1l}\\ ⋮& ⋰& ⋮& ⋰& & & & & \ddots & ⋮\\ & & & & & & & & & \\ & & & & & & & & & \\ ⋮& \ddots & ⋮& \ddots & & & & & ⋰& ⋮\\ a_{n1,k1}& \cdots & a_{n2,k1}& \cdots & & & & & \cdots & a_{nm,\mathrm{kl}}\end{array}\right)}.\end{array}$

[0110] However, this is a distorted way of viewing the problem because the internal coefficients a_νμ,στ are functioning on a deeper level than the external coefficients a_νμ. In more concrete terms, as mentioned above the solution to the linear prediction problem corresponding to aε(n,m,A) has nothing whatsoever to do with the linear prediction problem corresponding to aε(nk,ml,B).

[0111] The correct metaphor is to regard the expression (n,m,−) as defining a matrix class in the sense of object-oriented programming, then for any object A, (n,m,A) is an object inheriting the properties of (n,m,−), and utilizing the arithmetic of A to define operations such as matrix multiplication and addition. A itself inherits from a general scalar class defining the arithmetic of A. However, these classes are so general that (n,m,A) itself can be regarded as a scalar object, using its defined arithmetic. Accordingly, in the other direction, the scalar object A might itself be some matrix object (k,l,B).

[0112] In spite of the degree of abstraction this metaphor requires, it is the only one which correctly captures the general multi-channel situation. It is easy to imagine real-world multi-channel situations, such as the geophysics situation described previously, in which deep inheritance hierarchies are generated.

[0113] The present invention, according to one embodiment, addresses special cases of this general data-structuring problem, in which the introduction of non-commutative algebra into signal processing is a major advance towards a solution of the general case. The reason that multi-channel linear prediction produces significant data compression is the large cross-channel and cross-time correlation. This implies a high degree of redundancy in the datasets which can be removed, thereby reducing the bandwidth requirements.

[0114] Correlations are introduced in mechanical finite-element systems by physical constraints of shape, boundary conditions, material properties, and the like as well as the inertia of components with mass. This is also true for animal/robotic motion whose strongest constraints are due the semi-rigid structure of bone or metal.

[0115] In fact, as noted previously, multi-channel data is actually steeped with correlations—which was not an issue for single-channel processing. For example, when a single-channel linear predictor has been able to reduce the prediction error of a signal to 0, this can be interpreted as a sign of highly successful compression: it is demonstrated that the channel is carrying a deterministic sum of damped exponentials whose values can be determined by locating the roots of the characteristic polynomial of the system. In reality, things are not this simple; in practice, one regards a “perfect” linear prediction as indicative of too many coefficients and reduces the model order accordingly. However, things are far more complicated for multi-channel analysis because a large number of “perfect” channels are used.

[0116] That part of ordinary calculus, of any number of real or complex variables, which goes beyond simple algebra, is based in the fact that is a metric space for which the compact sets are precisely the closed, bounded sets. The higher-dimensional spaces ⁿ,ⁿ inherit the same property. The algebra of , plus the simple geometric combinatorics of covering regions by boxes allow all of calculus, complex, analysis, Fourier series and integrals, and the rest to be built up in the standard manner from this compactness property of .

[0117] Topologically and metrically, the quaternion ring is simply ⁴; with careful use of quaternion algebra (especially the non-commutativity), the same development can be followed for . All the standard results such as the Cauchy Integral Theorem, the Implicit Function Theorem, and the like have their quaternion analogs (often in left- and right-forms because of non-commutativity).

[0118] As a consequence, there is no problem in developing -versions of z-transforms and Laurent series, hence the P(z) and D(z) of the previous section. In fact, the theory of quaternion system functions is much richer than for the complex field because as is shown later, a quaternion variable z consists of two independent complex variables

[0119] $\left(\begin{array}{c}{z}_{+}\\ z_{-}\end{array}\right).$
Many unexpected frequency-domain phenomena will appear, unknown from the one variable situation, because of the geometric and analytic interactions of z₊ and z₋.

[0120] Because is non-commutative, the det( ) operator does not behave “properly”. The most important property of det( ) which fails over is its invariance under multiplication of columns or rows by a scalar; i.e., it is generally the case that

[0121] $\det \left(\begin{array}{ccccc}\begin{array}{c}{a}_{11}\\ ⋮\\ a_{M1}\end{array}& \cdots & k\left(\begin{array}{c}{a}_{1j}\\ a_{\mathrm{ij}}\\ a_{\mathrm{Mj}}\end{array}\right)& \cdots & \begin{array}{c}{a}_{1N}\\ a_{\mathrm{iN}}\\ a_{\mathrm{MN}}\end{array}\end{array}\right)\ne k\det \left(\begin{array}{ccccc}\begin{array}{c}{a}_{11}\\ ⋮\\ a_{M1}\end{array}& \cdots & \left(\begin{array}{c}{a}_{1j}\\ a_{\mathrm{ij}}\\ a_{\mathrm{Mj}}\end{array}\right)& \cdots & \begin{array}{c}{a}_{1N}\\ a_{\mathrm{iN}}\\ a_{\mathrm{MN}}\end{array}\end{array}\right),$
for kε.

[0122] As a result, basic identities such as det(ab)=det(a)det(b) and Cramer's Rule also fail.

[0123] Importantly, it is not the case that a matrix a over is invertible if and only if det(a) is invertible in . This is because the matrix adjoint a^adj generally satisfies a·a^adj≠det(a)·1 over non-commutative rings.

[0124] The present invention advantageously permits application of the Levinson algorithm in a wide class of cases in which the autocorrelation coefficients are not in a commutative field. In particular, it is shown that the modified Levinson algorithm applies to quaternion-valued autocorrelations, hence, for example, to 3 and (3+1)-dimensional data.

[0125] The algebra of complex numbers can be viewed as ordered pairs of real numbers (a,b), referred to as couplets. Addition was defined by the rule (a,b)+(c,d)=(a+c,b+d) and, most importantly, multiplication defined by the rule:
(a,b)·(c,d)=(ac−bd,ad−bc).

[0126] It has been shown that with these definitions, couplets could be added, subtracted, multiplied, and, when the divisor did not equal (0,0), divided as well.

[0127] Thus, i=√{square root over (−1)} can be simply defined as the couplet (0,1), while the couplet 1 (which is different in an abstract sense from the number 1) was defined to be (1,0).

[0128] Any couplet (a,b) could then be written uniquely in the form
(a,b)=a(1,0)+b(0,1)=a1+bi=a+bi
and the link to the complex numbers was complete.

[0129] An equivalent representation of the complex number a+bi is the (2×2) real matrix:

[0130] $〚a+\mathrm{bi}〛=\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right).$

[0131] This representation is important for understanding the more complicated quaternion representations.

[0132] Using the ordinary laws of matrix arithmetic, the following exists:

[0133] $〚a+\mathrm{bi}〛+〚c+\mathrm{di}〛=\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)+\left(\begin{array}{cc}c& d\\ -d& c\end{array}\right)=\left(\begin{array}{cc}a+c& b+d\\ -\left(b+d\right)& a+c\end{array}\right)=〚\left(a+\mathrm{bi}\right)+\left(c+\mathrm{di}\right)〛$ $\mathrm{and}$ $s·〚a+\mathrm{bi}〛=s·\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)=\left(\begin{array}{cc}s·a& s·b\\ -s·b& s·a\end{array}\right)=〚s·\left(a+\mathrm{bi}\right)〛,\mathrm{for}\mathrm{any}s\in R.$

[0134] Most significantly,

[0135] $\begin{array}{c}〚a+\mathrm{bi}〛·〚c+\mathrm{di}〛=\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)\left(\begin{array}{cc}c& d\\ -d& c\end{array}\right)\\ =\left(\begin{array}{cc}ac-\mathrm{bd}& \mathrm{ad}+\mathrm{bc}\\ -\left(\mathrm{ad}+\mathrm{bc}\right)& ac-\mathrm{bd}\end{array}\right)\\ =〚\left(a+\mathrm{bi}\right)·\left(c+\mathrm{di}\right)〛.\end{array}$

[0136] In this representation,

[0137] $1=〚1〛=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),I=〚i〛=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$
and thus

[0138] $\begin{array}{c}〚a+\mathrm{bi}〛=\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)=a·\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+b·\left(\begin{array}{cc}0& b\\ -b& 0\end{array}\right)=a·1+b·I\\ I^2=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)=\left(\begin{array}{cc}-1& 0\\ 0& -1\end{array}\right)=-1\end{array}$
and so, once again, the law i²=−1 receives a clear interpretation.

[0139] Also the complex conjugate is represented by the transpose:

[0140] $〚{\left(a+\mathrm{bi}\right)}^{*}〛=〚a-\mathrm{bi}〛=\left(\begin{array}{cc}a& -b\\ b& a\end{array}\right)={\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)}^T={〚a+\mathrm{bi}〛}^T$
and the squared norm |z|² represented by the determinant

[0141] ${a+\mathrm{bi}}^2=a^2+b^2=\det \left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)=\det 〚a+\mathrm{bi}〛.$

[0142] The following is noted:

[0143] $\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)·{\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)}^T=\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)·\left(\begin{array}{cc}a& -b\\ b& a\end{array}\right)=\left(a^2+b^2\right)·\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)=\left[\det \left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)\right]·\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$
and similarly

[0144] ${\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)}^T·\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)=\left[\det \left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)\right]·\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$

[0145] A real matrix C is called “orthogonal” if CC^T=C^TC=1, and the set of (n×n) real orthogonal matrices is denoted O(n). O(n) is a group under multiplication. A real matrix C is “extended orthogonal” if it satisfies the more general rule
CC^T=C^TC=r·1
for some rε and the set of (n×n) extended orthogonal matrices is denoted ⁺O(n). Thus, O(n)⊂⁺O(n). Since nr=trace(r·1)=trace(CC^T)≧0, where the trace of a matrix is the sum of the diagonal coefficients, r is necessarily non-negative and r=0C=0. So ⁺O(n)−{0} forms a group under matrix multiplication.

[0146] If C is orthogonal, then det(C)²=det(C)det(C^T)=det(CC^T)=det(1)=1 so det(C)=±1. An orthogonal matrix with det(C)=1 is called “special orthogonal,” and the set of (n×n) special orthogonal matrices (which is also a group) is denoted SO(n).

[0147] Analogously, an extended orthogonal matrix C is defined to be “special extended orthogonal” if det(C)≧0 and denote the set of special extended orthogonal matrices by S⁺O(n). Again SO(n)⊂S⁺O(n) and S⁺O(n)−{0} forms a group under multiplication.

[0148] It is observed that CεS⁺O(n) if and only if C=0 or (det(C)>0 and

[0149] $\frac{1}{{\hspace{0.17em}}^n\sqrt{\det \left(C\right)}}$
CεSO(n)). This implies that every CεS⁺O(n) has a unique representation C=sR, sε,s≧0, RεSO(n) and conversely. In particular, SO(n)={CεS⁺O(n)|det(C)=1}.

[0150] It can also be shown that a (2×2) real matrix C is special extended orthogonal if and only if it is of the form:

[0151] $C=\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right),a,b\in R,$
which are precisely the matrices with which represents. Thus this representation of is denoted by the S⁺O(2) representation.

[0152] In particular, the unit circle S¹={(x₁,x₂)ε²; x₁²+x₂²=1}≈{zε; |z|²=1} is isomorphic to the real rotation group SO(2) by means of the representation .

[0153] Instead of representing i by

[0154] $\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),$
it could be represented by

[0155] $\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right),$
and nothing in the arithmetic would differ. This is precisely the same phenomenon as in linear algebra in which it is more satisfactory in an abstract sense to define vector spaces merely by the laws they satisfy but in which computation is best performed in coordinate form by selecting some arbitrary basis.

[0156] A three-component analog of complex numbers (i.e., “triplets”) provides a useful arithmetic structure on three-dimensional space, just as the complex numbers put a useful arithmetic structure on two-dimensional space. The theory of addition and scalar multiplication for triplets, are as follows:
(a,b,c)+(d,e,ƒ)=(a+d,b+e,c+ƒ)
s·(a,b,c)=(s·a,s·b,s·c)

[0157] However, multiplying triplets is more difficult. Two ways of multiplication exist: dot product, cross product (i.e., vector product). The dot product (or the scalar product) is as follows:
(a,b,c)·(d,e,ƒ)=ad+be+cƒ
However, this product does not produce a triplet.

[0158] The other way is known as the cross product is as follows:
(a,b,c)×(d,e,ƒ)=(bƒ−ce,cd−aƒ,ae−bd).

[0159] The cross product has the advantage of producing a triplet from a pair of triplets, but fails to allow division. When A,B are triplets, the equation A×X=B is generally not solvable for X even when A≠0. However, the cross product contained the seed of the eventual solution in the anti-commutative law A×B=−B×A.

[0160] It is noted that three-dimensional space must be supplemented with a fourth temporal or scale dimension in order to form a complete system. Thus, 3-dimensional geometry must be embedded inside a (3+1)-dimensional geometry in order to have enough structure to allow certain types of objects (points at infinity, reciprocals of triplets, etc.) to exist.

[0161] The four-component objects named “quaternions,” have the usual addition and scalar multiplication laws. The definition of quaternion multiplication is as follows:
(a,b,c,d)·(e,ƒ,g,h)=(ae−bƒ−cg−dh,aƒ+be+ch−dg,ag+ce+dƒ−bh,ah+bg+de−cƒ)

[0162] Because of the complexity, this formula is not used for computation.

[0163] As with the representation of complex numbers as couplets, the first step is to define the units:
1=(1,0,0,0)
I=(0,1,0,0)
J=(0,0,1,0)
K=(0,0,0,1)

[0164] The previous formula then shows that I,J,K satisfy the multiplication rules:
I²=J²=K²=IJK=−1.

[0165] From these relations follow the permutation laws:
IJ=−JI=K
JK=−KJ=I
KI=−IK=J
and since 1a+Ib+Jc+Kd=(a,b,c,d)=a1+bI+cJ+cK, the usual laws of arithmetic combined with the above relations among the units defines quaternion multiplication completely. The quaternions is denoted as .

[0166] A quaternion has many representations, the most basic being the 4-vector form q=a1+bI+cJ+cK. Typically, the “1” is omitted (or identified with the number 1 where no ambiguity will result): q=a+bI+cJ+cK.

[0167] q=a+bI+cJ+cK naturally decomposes into its scalar part Sc(q)=aε and its vector (or principal) part Vc(q)=(bI+cJ+dK)ε³, where the quaternion units I,J,K are regarded as unit vectors in ³ forming a right-hand orthogonal basis.

[0168] q=Sc(q)+Vc(q) always holds. The expression, q=a+{right arrow over (ν)}, is used to indicate Sc(q)=a and Vc(q)={right arrow over (ν)}. This can be referred to as the (3+1)-vector representation of a quaternion.

[0169] The addition and scalar multiplication laws in the (3+1) form are simply
(a+{right arrow over (ν)})+(b+{right arrow over (w)})=(a+b)+({right arrow over (ν)}+{right arrow over (w)})
s·(a+{right arrow over (ν)})=(s·a+s·{right arrow over (ν)}),sε

[0170] However, the quaternion multiplication law in (3+1) form reveals the deep connection to the structure of three-dimensional space:
(a+{right arrow over (ν)})·(b+{right arrow over (w)})=(ab−{right arrow over (ν)}•{right arrow over (w)})+(a{right arrow over (w)}+b{right arrow over (ν)})+({right arrow over (ν)}×{right arrow over (w)}).

[0171] In the above expression, {right arrow over (ν)}•{right arrow over (w)} denotes dot product (cI+dJ+eK)═(ƒI+gJ+hK)=(cƒ+dg+eh) while {right arrow over (ν)}×{right arrow over (w)} denotes cross product

[0172] $\left(\mathrm{cI}+\mathrm{dJ}+\mathrm{eK}\right)\times \left(\mathrm{fI}+\mathrm{gJ}+\mathrm{hK}\right)=\begin{array}{ccc}c& f& I\\ d& g& J\\ e& h& K\end{array}=\left(\mathrm{dh}-\mathrm{eg}\right)I+\left(\mathrm{ef}-\mathrm{ch}\right)J+\left(\mathrm{cg}-\mathrm{df}\right)K.$

[0173] Since ab is ordinary scalar multiplication and a{right arrow over (w)},b{right arrow over (ν)} are just ordinary multiplications of a vector by a scalar, it can be seen that quaternion multiplication contains within it all four ways in which a pair of (3+1)-vectors can be multiplied.

[0174] It is suggestive that if the two relativistic spacetime intervals (Δx₁,Δy₁,Δz₁,cΔt₁),(Δx₂,Δy₂,Δz₂,cΔt₂) is represented by the quaternions
Δq₁=cΔt₁+(Δx₁)I+(Δy₁)J+(Δz₁)K,
Δq₂=cΔt₂+(Δx₂)I+(Δy₂)J+(Δz₂)K
then
Sc(Δq₁·Δq₂)=c²(Δt₁Δt₂)−(Δx₁Δx₂+Δy₁Δy₂+Δz₁Δz₂),
the familiar Minkowski scalar product.

[0175] The (3+1) product formula also shows that for any pure vector {right arrow over (ν)}, {right arrow over (ν)}²=−|{right arrow over (ν)}|²ε. In particular, when {circumflex over (ν)} is an ordinary unit vector in 3-space, {circumflex over (ν)}²=−1, which generalizes the rules for I,J,K.

[0176] As with the complex numbers, quaternions have a conjugation operation q*:
q*=(a+bI+cJ+dK)*=(a−bI−cJ−dK).

[0177] In (3+1) form this is (a+{right arrow over (ν)})*=(a−{right arrow over (ν)}). Generalizing the -formulae

[0178] ${\left(z^{*}\right)}^{*}=z,\mathrm{Re}\left(z\right)=\frac{1}{2}\left(z+z^{*}\right),i\mathrm{Im}\left(z\right)=\frac{1}{2}\left(z-z^{*}\right)$
yields the following:
(q*)*=q

[0179] $\begin{array}{c}Sc\left(q\right)=\frac{1}{2}\left(q+q^{*}\right).\\ Vc\left(q\right)=\frac{1}{2}\left(q-q^{*}\right)\end{array}$

[0180] Quaternions also have a norm generalizing the complex |z|=√{square root over (zz*)}:
|q|=√{square root over (qq*)}=√{square root over (q*q)}=√{square root over ((a²+b²+c²+d²))}ε
and, as with , |q|²≧0 and (|q|=0q=0). In (3+1) form the norm is calculated by |a+{right arrow over (ν)}|=√{square root over (a²+{right arrow over (ν)}•{right arrow over (ν)})}.

[0181] A unit quaternion is defined to be a uε such that |u|=1. It is noted that the quaternion units ±1,±I,±J,±K are all unit quaternions.

[0182] The chief peculiarity of quaternion arithmetic is the failure of the commutative law: for quaternions q,r, whereby generally q·r≠r·q; even the units do not commute: I·J=−J·I, etc. The (3+1) form (a+{right arrow over (ν)})·(b+{right arrow over (w)})=(ab−{right arrow over (ν)}•{right arrow over (w)})+(a{right arrow over (w)}+b{right arrow over (ν)})+({right arrow over (ν)}×{right arrow over (w)}) shows this most clearly. All the multiplication operations in this expression are commutative except the cross product {right arrow over (ν)}×{right arrow over (w)} which satisfies {right arrow over (ν)}×{right arrow over (w)}=−{right arrow over (w)}×{right arrow over (ν)}, hence is the source of non-commutativity. This also shows that if Vc(q) and Vc(r) are parallel vectors in ³ then q·r=r·q.

[0183] An important formula is the anti-commutative conjugate law
(q·r)*=r*·q*
which is most easily proved in the (3+1) form. Combined with the previous law (q*)*=q, this shows that conjugation is an anti-involution of .

[0184] Recall that the reciprocal of a non-zero complex number z can be written in the form

[0185] $z^{-1}=\frac{z^{*}}{{z}^2}$
and this also holds for quaternions:

[0186] $q^{-1}=\frac{q^{*}}{{q}^2},q\ne 0$
as is apparent by the calculation

[0187] $q\left(\frac{q^{*}}{{q}^2}\right)=\frac{q{q}^{*}}{{q}^2}=\frac{{q}^2}{{q}^2}=1$
and similarly for

[0188] $\left(\frac{q^{*}}{{q}^2}\right)q.$

[0189] As with all non-commutative groups, inverses anti-commute
(q≠0,r≠0)((qr)⁻¹=r⁻¹q⁻¹).

[0190] So possesses the four basic arithmetic operations but has a non-commutative multiplication, which is the definition of what is called a division ring.

[0191] A known result of Frobenius states that the only division rings which are finite-dimensional extensions of are itself (one-dimensional), the complex numbers (two-dimensional), and the quaternions ((3+1)-dimensional). This is another example of the exceptional properties of (3+1)-dimensional space.

[0192] The (n×n) identity matrix

[0193] $\hspace{1em}\left(\begin{array}{ccccc}1& 0& \cdots & & 0\\ 0& ⋰& & & ⋮\\ ⋮& & 1& & ⋮\\ ⋮& & & ⋰& 0\\ 0& & \cdots & 0& 1\end{array}\right)$
is denoted 1 to avoid confusion with the quaternion unit I.

[0194] There are many notations for the quaternion units; e.g., i,j,k; î,ĵ,{circumflex over (k)}; and I,J,K. A more general definition of the quaternions, based on is obtained as follows:

[0195] Let be a commutative field and e,ƒ,gε−{0}. (,e,ƒ,g), the quaternions over , is defined as the smallest -algebra which contains elements I,J,Kε(,e,ƒ,g) satisfying the relations
I²=−eƒ,J²=−eg,K²=−ƒg,IJK=−eƒg.

[0196] It can then be shown that
IJ=−JI=eK
JK=−KJ=gI.
KI=−IK=ƒJ

[0197] Any qε(,e,ƒ,g) can be written uniquely in the form q=a+bI+cJ+dK, a,b,c,dε with conjugate q*=a−bI−cJ−dK and norm ²|q|=a²+eƒb²+egc²+ƒgd².

[0198] An interesting situation is when the quadratic form w²+eƒx²+egy²+ƒgz² over is definite; i.e., (w²+eƒx²+egy²+ƒgz²=0)(w=x=y=z=0). In particular, for this to hold, none of −eƒ,−eg,−ƒg can be squares in . In this case, (,e,ƒ,g) is a division ring as well as a four-dimensional -algebra.

[0199] (,1,1,1)= are just Hamilton's quaternions.

[0200] In order to show that (,e,ƒ,g) exists, it is noted that the typical polynomial algebra constructions fail because the non-commutativity of the quaternion units.

[0201] Let be a -algebra, then the tensor algebra of over is the graded -algebra

[0202] $T_k\left(A\right)=\coprod _{n\ge 0}{\left(A{\otimes }_k\dots {\otimes }_{k}A\right)}_{n\mathrm{factors}}$
with product defined on basis elements by
(a₁{circumflex over (x)} . . . {circle around (x)}a_m)×(b₁{circle around (x)} . . . {circle around (x)} b_n)=(a₁{circle around (x)} . . . {circle around (x)}a_m{circle around (x)}b₁{circle around (x)} . . . {circle around (x)}b_n).

[0203] It is noted ({circle around (x)}_k . . . {circle around (x)}_k)_{0 factors}= by definition.

[0204] For e,ƒ,gε−{0}, define the quaternion -algebra (,e,ƒ,g) to be

[0205] $ℍ\left(𝕜,e,f,g\right)=\frac{T_{𝕜}\left({𝕜}^3\right)}{\Theta \left(𝕜,e,f,g\right)},$
where, defining I=(1,0,0),J=(0,1,0),K=(0,0,1), Θ(,e,ƒ,g) is the two-sided ideal generated by
eƒ+{circle around (x)}I
eg+J{circle around (x)}J
ƒg+K{circle around (x)}K
eƒg+I{circle around (x)}J{circle around (x)}K

[0206] The quaternion units {±1,±I,±J,±K} form a non-abelian group of order 8 under multiplication. By expressing as {1,1′,I,I′,J,J′,K,K′}, then the quaternions over any commutative field can be abstractly represented as the quotient ()=[]/Θ, where [] is the group ring and Θ is the two-sided ideal generated by 1+1′,I+I′,J+J′,K+K′.

[0207] There are many extensions ⊃ which are fields. For example, the field of formal quotients

[0208] $\frac{a_0+a_{1}x+\dots +a_n{x}^n}{b_0+b_{1}x+\dots +b_m{x}^m},$
a₀, a₁, . . . , a_n, b₀, b₁, . . . , b_mε. However, Frobenius' Theorem asserts that none of these can be finite-dimensional as vector spaces over .

[0209] Just as there are S⁺O(2) representations for the complex numbers, there are comparable representations for the quaternions. These are especially important because there are certain procedures, such as extracting the eigenstructure of quaternion matrices, which are nearly impossible except in these representations.

[0210] It is noted that an (n×n) complex matrix Q is called unitary if QQ*=Q*Q=1. Q* denotes the conjugate transpose also called the hermitian conjugate (which is sometimes denoted Q^H):

[0211] ${\left(\begin{array}{ccccc}{z}_{11}& & \cdots & & z_{1n}\\ & ⋰& & \ddots & \\ ⋮& & z_{ij}& & ⋮\\ & \ddots & & ⋰& \\ z_{\mathrm{n1}}& & \cdots & & z_{nn}\end{array}\right)}^{*}=\left(\begin{array}{ccccc}{z}_{11}^{*}& & \cdots & & z_{\mathrm{n1}}^{*}\\ & ⋰& & \ddots & \\ ⋮& & z_{ij}^{*}& & ⋮\\ & \ddots & & ⋰& \\ z_{1n}^{*}& & \cdots & & z_{nn}^{*}\end{array}\right).$

[0212] It is noted when Q is real, Q*=Q^T. The group of (n×n) unitary matrices is denoted U(n). Thus O(n)⊂U(n).

[0213] As with the orthogonal matrices, a complex matrix Q is termed “extended unitary” if the more general rule
QQ*=Q*Q=r·1,rε
holds and denote the (n×n) extended unitary matrices by ⁺U(n). So ⁺O(n)∪⁺U(n)⊂⁺U(n) and ⁺U(n)−{0} is a group under multiplication.

[0214] A unitary matrix Q is special unitary if det(Q)=1 and analogously an extended unitary matrix Q is special extended unitary if det(Q)≧0. The special extended unitary matrices are denoted S⁺U(n); thus, (S⁺O(n)∪SU(n))⊂S⁺U(n), and S⁺U(n)−{0} is a group under multiplication.

[0215] As with S⁺O(n), it is straight forward to calculate that QεS⁺U(n) if and only if Q=0 or (det(Q)ε,det(Q)>0 and

[0216] $\frac{1}{\sqrt[n]{\det \left(Q\right)}}$
QεSU(n)). This implies that every QεS⁺U(n) has a unique representation Q=sU, sε,s≧0, UεSU(n) and conversely.

[0217] It can be shown that a (2×2) complex matrix Q is special extended unitary if and only if it is of the form:

[0218] $Q=\left(\begin{array}{cc}{z}_{+}& z_{-}\\ -z_{-}^{*}& z_{+}^{*}\end{array}\right),z_{+},z_{-}\in C.$

[0219] Defining

[0220] $〚z_{+}+z_{-}J〛=\left(\begin{array}{cc}{z}_{+}& z_{-}\\ -z_{-}^{*}& z_{+}^{*}\end{array}\right),$
it can be shown, using the laws of quaternion arithmetic in the bicomplex representation, that converts all the algebraic operations in into matrix operations. is called the S⁺U(2) representation.

[0221] Moreover, the S⁺U(2) representation sends conjugation to hermitian conjugation and the squared norm to the determinant:

[0222] $〚{\left(z_{+}+z_{-}J\right)}^{*}〛=〚z_{+}^{*}-z_{-}J〛=\left(\begin{array}{cc}{z}_{+}^{*}& -z_{-}\\ z_{-}^{*}& z_{+}\end{array}\right)={\left(\begin{array}{cc}{z}_{+}& z_{-}\\ -z_{-}^{*}& z_{+}^{*}\end{array}\right)}^{*}={〚z_{+}+z_{-}J〛}^{*}z_{+}+z_{-}J{}^2=z_{+}{}^2+z_{-}{}^2=\det \left(\begin{array}{cc}{z}_{+}& z_{-}\\ -z_{-}^{*}& z_{+}^{*}\end{array}\right)=\det 〚z_{+}+z_{-}J〛.$
In particular, the unit 3-sphere
S³={(x₁,x₂,x₃,x₄)ε⁴;x₁²+x₂²+x₃²+x₄²=1}≈{qε;|q|²=1}
is isomorphic to the spin group SU(2) by means of the representation .

[0223] The unit quaternions {qε; |q|²=1} is denoted ⊂. In terms of the (3+1) form of quaternions, the S⁺U(2) representation is

[0224] $〚a+\mathrm{bI}+\mathrm{cJ}+\mathrm{cK}〛=\left(\begin{array}{cc}a+\mathrm{bi}& c+\mathrm{di}\\ -c+\mathrm{di}& a-\mathrm{bi}\end{array}\right).$

[0225] Decomposing the matrix a+bI+cJ+cK yields

[0226] $\begin{array}{c}〚a+\mathrm{bI}+\mathrm{cJ}+\mathrm{cK}〛=\left(\begin{array}{cc}a+\mathrm{bi}& c+\mathrm{di}\\ -c+\mathrm{di}& a-\mathrm{bi}\end{array}\right)\\ =a\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+b\left(\begin{array}{cc}i& 0\\ 0& -i\end{array}\right)+c\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)+d\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)\end{array}$
and, thus,

[0227] $〚1〛=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),〚I〛=\left(\begin{array}{cc}i& 0\\ 0& -i\end{array}\right),\hspace{1em}〚J〛=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),〚K〛=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right).$

[0228] The above are denoted as the standard units of the S⁺U(2) representation.

[0229] It is also easy to extend the S⁺U(2) representation to m×n quaternion matrices component wise:

[0230] $〚\left(\begin{array}{cccc}& ⋮& & \\ \dots & a_{\mathrm{ij}}& \dots & \dots \\ & ⋮& & \end{array}\right)〛=\left(\begin{array}{cccc}& ⋮& & \\ \dots & 〚a_{\mathrm{ij}}〛& \dots & \dots \\ & ⋮& & \end{array}\right).$

[0231] This representation will preserve all the additive and multiplicative properties of quaternion matrices.

[0232] Assuming {circumflex over (α)}ε³ is a unit vector and θε be an angle, then the quaternion is defined as follows:

[0233] $u=u\left(\theta ,\hat{\alpha }\right)=\cos \frac{\theta }{2}+\left(\sin \frac{\theta }{2}\right)\hat{\alpha }.$

[0234] For all vectors {right arrow over (ν)}ε³, the quaternion product u{right arrow over (ν)}u* is also a vector and is the right-handed rotation of {right arrow over (ν)} about the axis {circumflex over (α)} by angle θ. It is noted u(θ,{circumflex over (α)}) is always a unit quaternion; i.e., u(θ,{circumflex over (α)})ε.

[0235] This result has found uses in, for example, computer animation and orbital mechanics because it reduces the work required to compound rotations: a series of rotations (θ₁,{circumflex over (α)}₁), . . . ,(θ_k,{circumflex over (α)}_k) can be represented by the quaternion product u(θ_k,{circumflex over (α)}_k) . . . u(θ₁,{circumflex over (α)}₁) which is much more efficient to compute than the product of the associated rotation matrices. Moreover, by inverting the map (θ,{circumflex over (α)})(θ,{circumflex over (α)}) the resultant angle and axis of this series of rotations can be calculated:
(θ_net,{circumflex over (α)}_net)=u⁻¹[u(θ_k,{circumflex over (α)}_k) . . . u(θ₁,{circumflex over (α)}₁)],
which is simpler than computing the eigenstructure of the product rotation matrix.

[0236] If q=a+{right arrow over (ν)} is an arbitrary quaternion and uε then uqu*=u(a+{right arrow over (ν)})u*=auu*+u{right arrow over (ν)}u*=a+u{right arrow over (ν)}u* so that rotation by u leaves Sc(q) unchanged. In particular, when qε, uqu*=q so rotation leaves ⊂ invariant. Thus ulu*=1.

[0237] Also
u(q+r)u*=uqu*+uru*
u(qr)u*=u(q(u*u)r)u*=(uqu*)(uru*)
(uqu*=r)(q=u*ru).

[0238] The conclusion is that the rotation map q(uqu*) is an algebraic automorphism of i.e., a structure-preserving one-to-one correspondence.

[0239] Assuming {right arrow over (u)},{right arrow over (ν)} are non-parallel vectors of the same length, then there is at least one rotation of ³ which sends {right arrow over (u)} to {right arrow over (ν)}. Any unit vector {circumflex over (α)} which lies on the plane of points which are equidistant from the tips of {right arrow over (u)},{right arrow over (ν)} can be used as an axis for a rotation which sends {right arrow over (u)} to {right arrow over (ν)}.

[0240] As {right arrow over (u)} is rotated around one of these axes, the tip of {right arrow over (u)} moves in a circle which lies in the sphere centered at the origin and passing through the tips of {right arrow over (u)},{right arrow over (ν)}. Generally this is a small circle on this sphere. However, there are two unit vectors {circumflex over (α)} around which the tip of {right arrow over (u)} moves in a great circle; namely

[0241] $\hat{\alpha }=\pm \frac{\overrightarrow{u}\times \overrightarrow{v}}{\overrightarrow{u}\times \overrightarrow{v}},$
the unique unit vectors perpendicular to both {right arrow over (u)} and {right arrow over (ν)}.

[0242] When rotated around such {circumflex over (α)}, the tip of {right arrow over (u)} moves along either the longest or shortest path between the tips depending on the orientations. In either case, this path is an external of the length of the paths. Any unit vector around which {right arrow over (u)} can be rotated into {right arrow over (ν)} along an external path is referred to as an “external unit vector.” Clearly {circumflex over (α)} is an external unit vector, then so is −{circumflex over (α)}.

[0243] When {right arrow over (u)}={right arrow over (ν)}≠{right arrow over (0)}, the external vectors are

[0244] $\hat{\alpha }=\pm \frac{\overrightarrow{u}}{\overrightarrow{u}}$
since any rotation fixing {right arrow over (u)} must have the line containing {right arrow over (u)} as an axis. When {right arrow over (u)}=−{right arrow over (ν)}≠{right arrow over (0)}, the external vectors are all unit vectors in the plane perpendicular to {right arrow over (u)}. When {right arrow over (u)}={right arrow over (ν)}={right arrow over (0)}, the external vectors are all unit vectors.

[0245] Now, it is assumed that {circumflex over (α)},{circumflex over (β)},{circumflex over (γ)} and {circumflex over (α)}′,{circumflex over (β)}′,{circumflex over (γ)}′ are two right-handed, orthonormal systems of vectors: {circumflex over (α)}⊥{circumflex over (β)}, |{circumflex over (α)}|={circumflex over (β)}|=1, {circumflex over (γ)}={circumflex over (α)}×{circumflex over (β)} and similarly for {circumflex over (α)}′,{circumflex over (β)}′,{circumflex over (γ)}′. To simplify the analysis, that it is further assumed that {circumflex over (α)},{circumflex over (α)}′ are not parallel and {circumflex over (β)},{circumflex over (β)}′ are not parallel.

[0246] As discussed above, all the rotations sending {circumflex over (α)} to {circumflex over (α)}′ determine a plane and similarly for the rotations sending {circumflex over (β)} to {circumflex over (β)}′. Assuming these planes are not the same, they will intersect in a line through the origin. There is then a unique rotation around this line (and only around this line) which will simultaneously send {circumflex over (α)} to {circumflex over (α)}′ and {circumflex over (β)} to {circumflex over (β)}′. Since {circumflex over (γ)}={circumflex over (α)}×{circumflex over (β)} and {circumflex over (γ)}′={circumflex over (α)}′×{circumflex over (α)}′, this rotation also sends {circumflex over (γ)} to {circumflex over (γ)}′.

[0247] By carefully analyzing the various cases when parallelism occurs, the following can be shown:

[0248] Proposition 1 For any two right-handed, orthonormal systems of vectors {circumflex over (α)},{circumflex over (β)},{circumflex over (γ)} and {circumflex over (α)}′,{circumflex over (β)}′,{circumflex over (γ)}′, there is a unit quaternion uε such that
{circumflex over (α)}′=u{circumflex over (α)}u*
{circumflex over (β)}′=u{circumflex over (β)}u*,
{circumflex over (γ)}′=u{circumflex over (α)}u*
Moreover, u is unique up to sign: ±u will both work.

[0249] The sign ambiguity is easy to understand:

[0250] $u=u\left(\theta ,\hat{\alpha }\right)=\cos \frac{\theta }{2}+\left(\sin \frac{\theta }{2}\right)$
{circumflex over (α)} is the rotation around {circumflex over (α)} by angle θ while

[0251] $\begin{array}{c}-u=-\cos \frac{\theta }{2}-\left(\sin \frac{\theta }{2}\right)\hat{\alpha }\\ =\cos \left(\frac{2\pi -\theta }{2}\right)+\sin \left(\frac{2\pi -\theta }{2}\right)\left(-\hat{\alpha }\right)\\ =u\left(\left(2\pi -\theta \right),-\hat{\alpha }\right)\end{array}$
is the rotation around −{circumflex over (α)} by angle (2π−θ). However, these are geometrically identical operations.

[0252] Because of the automorphism properties, if uε and the following is defined
I′=uIu*
J′=uJu*
K′=uKu*
then the relations
I′²=J′²=K′²=I′J′K′=−1
I′J′=K′,J′K′=I′,K′I′=J′
will hold. This means the new units I′,J′,K′ are algebraically indistinguishable form the old units I,J,K.

[0253] Therefore, any right-handed, orthonormal system of unit vectors can function as the quaternion units.

[0254] As a result of this, neither the bicomplex nor the S⁺U(2) representations are unique. For example, it was mentioned previously that any of the maps
(a+bi)(a+bI)
(a+bi)(a+bJ)
(a+bi)(a+bK)
could be used to define a distinct embedding ⊂ hence induces a distinct bicomplex representation of .

[0255] All of these arise by cyclically permuting the units: I,J,K→J,K,I→K,I,J which can be accomplished by the rotation quaternion

[0256] $u=\frac{1}{\sqrt{3}}\left(I+J+K\right).$
In fact, there are exactly 24 different right-hand systems that can be selected from {±I,±J,±K}, any of which can function as a quaternion basis, and all of which are obtained by some rotation quaternion of the form

[0257] $u=\frac{1}{\sqrt{3}}\left(I\pm J\pm K\right).$

[0258] In other words, if UεSU(2), then

[0259] $\begin{array}{c}{〚a+bI+cJ+cK〛}_U=a\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+b\left[U\left(\begin{array}{cc}i& 0\\ 0& -i\end{array}\right)U^{*}\right]+\\ c\left[U\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)U^{*}\right]+\\ d\left(U\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)U^{*}\right]\end{array}$
is a valid S⁺U(2) representation.

[0260] This illustrates the additional richness of the quaternions over the complex numbers: the only non-trivial -invariant automorphism of is complex conjugation but has a distinct automorphism for each unit {±u}⊂.⁶

[0261] Assuming a is an n×n matrix over . a is called normal if it commutes with its conjugate: aa*=a*a. Important classes of normal matrices include the following:

[0262] Hermitian (or symmetric or self-adjoint): a*=a

[0263] Anti-hermitian (or anti-symmetric): a*=−a

[0264] Unitary (or orthogonal): a*=a⁻¹

[0265] Non-negative: a=bb* for some b

[0266] Semi-positive: a is non-negative and a≠0

[0267] A projection: a²=a*=a

[0268] It is a classic result that any normal matrix a can be diagonalized by a unitary matrix; that is, there is a unitary matrix u and a diagonal matrix

[0269] $\lambda =\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & ⋰& \\ & & & {\lambda }_n\end{array}\right)$
such that u*au=λ.

[0270] λ₁, λ₂, . . . , λ_nε are the eigenvalues of a and the columns of u form an orthonormal basis for ⁿ with the inner product

[0271] $\langle \overrightarrow{x},\overrightarrow{y}\rangle =\sum _v{x}_v{y}_v^{*}.$

[0272] The standard normal classes can be characterized by the properties of λ₁, λ₂, . . . , λ_n:

[0273] Hermitianλ₁, λ₂, . . . , λ_nε

[0274] Anti-hermitian

[0275] $\frac{1}{i}{\lambda }_1,\frac{1}{i}{\lambda }_2,\dots ,\frac{1}{i}{\lambda }_n,\in R$

[0276] Unitary|λ₁|=|λ₂|= . . . =|λ_n|=1

[0277] Non-negativeλ₁, λ₂, . . . , λ_nε and λ₁, λ₂, . . . , λ_n≧0

[0278] Semi-positiveλ₁, λ₂, . . . , λ_nε and for some ν, λ_ν>0.

[0279] A projectionλ₁, λ₂, . . . , λ_nε{0,1}

[0280] In particular, it is noted that any real normal matrix aε^n×n will generally have complex eigenvalues and eigenvectors. In the special case that a is symmetric (a^T=a), a can be diagonalized by a real orthogonal matrix and has real diagonal entries.

[0281] The first step in quaternion modeling is to generalize this result to ; i.e., to show that any normal quaternion matrix a can be diagonalized by a unitary quaternion matrix. In fact, it can be shown that the eigenvalues are in ⊂. This latter fact is important because it means the characteristic polynomial p_a(λ)=det(λ1−a) need not be discussed, which, as mentioned above, is badly behaved over . This also implies that the same classification of the normal types based on the properties of λ₁, λ₂, . . . , λ_nε works for quaternion matrices as well.

[0282] This can be regarded as the Fundamental Theorem of quaternions because it has so many important consequences. In particular, in the case n=1, this will yield the polar representation of a quaternion, which is the basis for quaternion spatial modeling.

[0283] As pointed out above, parts of standard linear algebra do not work over . However, linear independence and the properties of span( ) in ⁿ work the same way as in ⁿ except that the left scalar multiplication needs to be distinguished from the right scalar multiplication. Because is a division ring, the following lemmas result:

[0284] Lemma 1 Let {right arrow over (w)}, {right arrow over (ν)}₁, . . . , {right arrow over (ν)}_lεⁿ and suppose {{right arrow over (ν)}₁, . . . , {right arrow over (ν)}_l} is linearly independent but {{right arrow over (w)}, {right arrow over (ν)}₁, . . . , {right arrow over (ν)}_l} is linearly dependent, then {right arrow over (w)}_kεspan({right arrow over (ν)}₁, . . . , {right arrow over (ν)}_l).

[0285] Lemma 2 Let {right arrow over (w)}₁, . . . , {right arrow over (w)}_k, {right arrow over (ν)}₁, . . . , {right arrow over (ν)}_lεⁿ such that {right arrow over (w)}₁, . . . , {right arrow over (w)}_kεspan({right arrow over (ν)}₁, . . . , {right arrow over (ν)}_l) and k>l, then {{right arrow over (w)}₁, . . . , {right arrow over (w)}_k} is linearly dependent.

[0286] These lemmas imply all the usual results concerning bases and dimension including the fact that any linearly independent set can be extended to a basis for ⁿ.

[0287] The inner product yields:

[0288] $\langle \overrightarrow{x},\overrightarrow{y}\rangle =\langle \left(\begin{array}{c}\begin{array}{c}{x}_1\\ ⋮\end{array}\\ x_n\end{array}\right),\left(\begin{array}{c}\begin{array}{c}{y}_1\\ ⋮\end{array}\\ y_n\end{array}\right)\rangle =\underset{v=1}{\sum ^n}{x}_v{y}_v^{*}$
which satisfies the usual properties of the inner product over ⁿ including {right arrow over (x)},{right arrow over (x)}=0({right arrow over (x)}=0) and q{right arrow over (x)},{right arrow over (y)}=q·{right arrow over (x)},{right arrow over (y)}, qε. Perpendicularity is defined by ({right arrow over (x)}⊥{right arrow over (y)}){right arrow over (x)},{right arrow over (y)}=0.

[0289] Lemma 3 (Projection Theorem for ) Let {right arrow over (ν)}₁, . . . , {right arrow over (ν)}_lεⁿ, then for all {right arrow over (w)}εⁿ, there exist q₁, . . . q_lε and a unique {right arrow over (e)}εⁿ such that {right arrow over (w)}q₁{right arrow over (ν)}_l+ . . . q_l{right arrow over (ν)}_l+{right arrow over (e)} and {right arrow over (e)}⊥{right arrow over (ν)}₁, . . . , {right arrow over (ν)}_l. If {{right arrow over (ν)}₁, . . . , {right arrow over (ν)}_l} is linearly independent, then q₁, . . . q_l are also unique.

[0290] Using the Projection Theorem, it can be shown that ⁿ has an orthonormal basis and, in fact, any orthonormal set {{right arrow over (ν)}₁, . . . , {right arrow over (ν)}_l} can be extended to an orthonormal basis.

[0291] The matrix u of change-of-basis to any orthonormal set is unitary and thus the matrix g of any linear operator

[0292] $n^{}\stackrel{G}{->}{n}^{}$
is transformed to ugu* by the basis change.

[0293] Let

[0294] $a=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$
be a 2×2 matrix over . Define the matrix

[0295] $a^{\diamond }=\left(\begin{array}{cc}{d}^{*}& -c^{*}\\ -b^{*}& a^{*}\end{array}\right)$
and suppose

[0296] $a\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{c}x\\ y\end{array}\right),\mathrm{then}{a}^{\diamond }\left(\begin{array}{c}-v^{*}\\ u^{*}\end{array}\right)=\left(\begin{array}{cc}{d}^{*}& -c^{*}\\ -b^{*}& a^{*}\end{array}\right)\left(\begin{array}{c}-v^{*}\\ u^{*}\end{array}\right)=\left(\begin{array}{c}-{\left(\mathrm{cu}+\mathrm{dv}\right)}^{*}\\ {\left(au+\mathrm{bv}\right)}^{*}\end{array}\right)=\left(\begin{array}{c}-y^{*}\\ x^{*}\end{array}\right).$
Next it is noted that for any

[0297] $\left(\begin{array}{cc}{z}_{+}& z_{-}\\ -z_{-}^{*}& z_{+}^{*}\end{array}\right)\in S^{+}U\left(2\right),{\left(\begin{array}{cc}{z}_{+}& z_{-}\\ -z_{-}^{*}& z_{+}^{*}\end{array}\right)}^0=\left(\begin{array}{cc}{\left(z_{+}^{*}\right)}^{*}& -{\left(-z_{-}^{*}\right)}^{*}\\ -{\left(z_{-}\right)}^{*}& {\left(z_{+}\right)}^{*}\end{array}\right)=\left(\begin{array}{cc}{z}_{+}& z_{-}\\ -z_{-}^{*}& z_{+}^{*}\end{array}\right).$
Thus, the following lemma results:

[0298] Lemma 4 Let qε and

[0299] $\left(\begin{array}{c}u\\ v\end{array}\right),\left(\begin{array}{c}x\\ y\end{array}\right)\in C^2$
such that

[0300] $〚q〛\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{c}x\\ y\end{array}\right),$
then

[0301] $〚q〛\left(\begin{array}{c}-v^{*}\\ u^{*}\end{array}\right)=\left(\begin{array}{c}-y^{*}\\ x^{*}\end{array}\right).$

[0302] It is noted that this result is independent of which form of is used. However, the next result requires selecting a specific form:

[0303] Proposition 2 It is assumed that a be an n×n quaternion matrix and {right arrow over (w)}ε²ⁿ−{{right arrow over (0)}} is an eigenvector of the standard representation a with eigenvalue λε, {right arrow over (w)} can be written in the form

[0304] $\overrightarrow{w}=\left(\begin{array}{c}{u}_1\\ v_1\\ ⋮\\ u_n\\ v_n\end{array}\right).$
Also, λε can be identified with λε by replacing iε by Iε; then

[0305] $a\left(\begin{array}{c}{u}_1-{\mathrm{Jv}}_1\\ ⋮\\ u_n-{\mathrm{Jv}}_n\end{array}\right)=\left(\begin{array}{c}{u}_1-{\mathrm{Jv}}_1\\ ⋮\\ u_n-{\mathrm{Jv}}_n\end{array}\right)·\lambda .$
Writing a and

[0306] $\overrightarrow{w}=\left(\begin{array}{c}{u}_1\\ v_1\\ ⋮\\ u_n\\ v_n\end{array}\right)$
in blocks as

[0307] $〚a〛=\left(\begin{array}{ccc}& & \\ & 〚a_{\mathrm{kl}}〛& \\ & & \end{array}\right)\mathrm{and}\stackrel{->}{w}=\left(\begin{array}{c}{u}_1\\ v_1\\ \\ \\ u_n\\ v_n\end{array}\right),$
the equation a{right arrow over (w)}={right arrow over (w)}λ is seen to be

[0308] $\sum _{l=1}^n〚a_{\mathrm{kl}}〛\left(\begin{array}{c}{u}_l\\ v_l\end{array}\right)=\left(\begin{array}{c}{u}_k\\ v_k\end{array}\right)\lambda =\left(\begin{array}{c}{u}_k\lambda \\ v_k\lambda \end{array}\right),$
k=1, . . . , n.

[0309] By Lem. 3,

[0310] $\sum _{l=1}^n〚a_{\mathrm{kl}}〛\left(\begin{array}{c}-v_l^{*}\\ u_l^{*}\end{array}\right)=\left(\begin{array}{c}-v_k^{*}{\lambda }^{*}\\ u_k^{*}{\lambda }^{*}\end{array}\right)=\left(\begin{array}{c}-v_k^{*}\\ u_k^{*}\end{array}\right)·{\lambda }^{*},k=1,\cdots ,n\Rightarrow \sum _{l=1}^n〚a_{\mathrm{kl}}〛\left(\begin{array}{cc}{u}_l& -v_l^{*}\\ v_l& u_l^{*}\end{array}\right)=\left(\begin{array}{cc}{u}_k& -v_k^{*}\\ v_k& u_k^{*}\end{array}\right)·\left(\begin{array}{cc}\lambda & 0\\ 0& {\lambda }^{*}\end{array}\right),k=1,\cdots ,n.\mathrm{However},\left(\begin{array}{cc}{u}_l& -v_l^{*}\\ v_l& u_l^{*}\end{array}\right)=\left(\begin{array}{cc}{u}_l& \left(-v_l^{*}\right)\\ -{\left(-v_l^{*}\right)}^{*}& u_l^{*}\end{array}\right)=u_l+\left(-v_l^{*}\right)J〛=〚u_l-J{v}_l〛\mathrm{and}\left(\begin{array}{cc}\lambda & 0\\ 0& {\lambda }^{*}\end{array}\right)=〚\lambda +0J〛=〚\lambda 〛$
in the standard representation.

[0311] Therefore

[0312] $\sum _{l=1}^n{a}_{\mathrm{kl}}\left(u_l-J{v}_l\right)=\left(u_k-J{v}_k\right)·\lambda \mathrm{in}$ $H\Rightarrow a\left(\begin{array}{c}{u}_1-J{v}_1\\ ⋮\\ u_n-J{v}_n\end{array}\right)=\left(\begin{array}{c}{u}_1-J{v}_1\\ ⋮\\ u_n-J{v}_n\end{array}\right)·\lambda \mathrm{in}{H}^{\underset{\_}{w}}.$

[0313] It is noted that this proposition shows that if column vectors are used to represent ⁿ then “eigenvalue” must be taken to mean “right eigenvalue”.

[0314] Proposition 3 (The Fundamental Theorem): Let a be an n×n normal matrix over , then there exists an n×n unitary matrix u over and a diagonal matrix

[0315] $\lambda =\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & ⋰& \\ & & & {\lambda }_n\end{array}\right)$
with λ₁, λ₂, . . . , λ_nε such that u*au=λ. λ is unique up to permutations of the diagonal coefficients.

[0316] Let a be normal. Since every matrix over ²ⁿ has an eigenvector, Prop. 2 implies that a has an eigenvector {right arrow over (y)}εⁿ−{{right arrow over (0)}} with eigenvalue λ₁ε. By the corollaries to the Projection Theorem, {right arrow over (γ)} can be extended to an orthogonal basis for ⁿ. In this basis, a becomes

[0317] $u_1^{*}a{u}_1=\left(\begin{array}{cccc}{\lambda }_1& q_2& \cdots & q_n\\ 0& & & \\ ⋮& & a^{\prime }& \\ 0& & & \end{array}\right),$
where u₁ is unitary. This matrix is also normal and since