k-Plane Clustering

A finite new algorithm is proposed for clustering m given points in n-dimensional real space into k clusters by generating k planes that constitute a local solution to the nonconvex problem of minimizing the sum of squares of the 2-norm distances between each point and a nearest plane. The key to the algorithm lies in a formulation that generates a plane in n-dimensional space that minimizes the sum of the squares of the 2-norm distances to each of m₁ given points in the space. The plane is generated by an eigenvector corresponding to a smallest eigenvalue of an n × n simple matrix derived from the m₁ points. The algorithm was tested on the publicly available Wisconsin Breast Prognosis Cancer database to generate well separated patient survival curves. In contrast, the k-mean algorithm did not generate such well-separated survival curves.     

QR factorization of Toeplitz matrices

Summary This paper presents a new algorithm for computing theQR factorization of anm×n Toeplitz matrix inO(mn) operations. The algorithm exploits the procedure for the rank-1 modification and the fact that both principal (m–1)×(n–1) submatrices of the Toeplitz matrix are identical. An efficient parallel implementation of the algorithm is possible.     

Open questions concerning Weiszfeld's algorithm for the Fermat-Weber location problem

The Fermat—Weber location problem is to find a point in ⁿ that minimizes the sum of the weighted Euclidean distances fromm given points in ⁿ . A popular iterative solution method for this problem was first introduced by Weiszfeld in 1937. In 1973 Kuhn claimed that if them given points are not collinear then for all but a denumerable number of starting points the sequence of iterates generated by Weiszfeld's scheme converges to the unique optimal solution. We demonstrate that Kuhn's convergence theorem is not always correct. We then conjecture that if this algorithm is initiated at the affine subspace spanned by them given points, the convergence is ensured for all but a denumerable number of starting points.     

Efficient Algorithms for the Smallest Enclosing Ball Problem

Consider the problem of computing the smallest enclosing ball of a set of m balls in ⁿ. Existing algorithms are known to be inefficient when n > 30. In this paper we develop two algorithms that are particularly suitable for problems where n is large. The first algorithm is based on log-exponential aggregation of the maximum function and reduces the problem into an unconstrained convex program. The second algorithm is based on a second-order cone programming formulation, with special structures taken into consideration. Our computational experiments show that both methods are efficient for large problems, with the product mn on the order of 10⁷. Using the first algorithm, we are able to solve problems with n = 100 and m = 512,000 in about 1 hour.His work was supported by Australian Research Council.Research supported in part by the Singapore-MIT Alliance.     

A summability method for the arithmetic Fourier transform

The Arithmetic Fourier Transform (AFT) is an algorithm for the computation of Fourier coefficients, which is suitable for parallel processing and in which there are no multiplications by complex exponentials. This is accomplished by the use of the Möbius function and Möbius inversion. However, the algorithm does require the evaluation of the function at an array of irregularly spaced points. In the case that the function has been sampled at regularly spaced points, interpolation is used at the intermediate points of the array. Generally theAFT is most effective when used to calculate the Fourier cosine coefficients of an even function.In this paper a summability method is used to derive a modification of theAFT algorithm. The proof of the modification is quite independent of theAFT itself and involves a summation by primes. One advantage of the new algorithm is that with a suitable sampling scheme low order Fourier coefficients may be calculated without interpolation.     

The Constrained Solutions of Two Matrix Equations

We study the symmetric positive semidefinite solution of the matrix equation AX ₁ A ^T + BX ₂ B ^T = C, where A is a given real m×n matrix, B is a given real m×p matrix, and C is a given real m×m matric, with m, n, p positive integers; and the bisymmetric positive semidefinite solution of the matrix equation D ^T XD = C, where D is a given real n×m matrix, C is a given real m×m matrix, with m, n positive integers. By making use of the generalized singular value decomposition, we derive general analytic formulae, and present necessary and sufficient conditions for guaranteeing the existence of these solutions. Received December 17, 1999, Revised January 10, 2001, Accepted March 5, 2001     

The Fermat—Weber location problem revisited

The Fermat—Weber location problem requires finding a point in ^N that minimizes the sum of weighted Euclidean distances tom given points. A one-point iterative method was first introduced by Weiszfeld in 1937 to solve this problem. Since then several research articles have been published on the method and generalizations thereof. Global convergence of Weiszfeld's algorithm was proven in a seminal paper by Kuhn in 1973. However, since them given points are singular points of the iteration functions, convergence is conditional on none of the iterates coinciding with one of the given points. In addressing this problem, Kuhn concluded that whenever them given points are not collinear, Weiszfeld's algorithm will converge to the unique optimal solution except for a denumerable set of starting points. As late as 1989, Chandrasekaran and Tamir demonstrated with counter-examples that convergence may not occur for continuous sets of starting points when the given points are contained in an affine subspace of ^N . We resolve this open question by proving that Weiszfeld's algorithm converges to the unique optimal solution for all but a denumerable set of starting points if, and only if, the convex hull of the given points is of dimensionN.     

Construction of Matrices with Prescribed Singular Values and Eigenvalues

Two issues concerning the construction of square matrices with prescribe singular values an eigenvalues are addressed. First, a necessary and sufficient condition for the existence of an n × n complex matrix with n given nonnegative numbers as singular values an m ( n) given complex numbers to be m of the eigenvalues is determined. This extends the classical result of Weyl and Horn treating the case when m = n. Second, an algorithm is given to generate a triangular matrix with prescribe singular values an eigenvalues. Unlike earlier algorithms, the eigenvalues can be arranged in any prescribe order on the diagonal. A slight modification of this algorithm allows one to construct a real matrix with specified real an complex conjugate eigenvalues an specified singular values. The construction is done by multiplication by diagonal unitary matrices, permutation matrices and rotation matrices. It is numerically stable and may be useful in developing test software for numerical linear algebra packages.     

Parallel Dixon Matrices by Bracket

It is known that the Dixon matrix can be constructed in parallel either by entry or by diagonal. This paper presents another parallel matrix construction, this time by bracket. The parallel by bracket algorithm is the fastest among the three, but not surprisingly it requires the highest number of processors. The method also shows analytically that the Dixon matrix has a total of m(m+1)²(m+2)n(n+1)²(n+2)/36 brackets but only mn(m+1)(n+1)(mn+2m+2n+1)/6 of them are distinct.     

Approximation algorithms for scheduling unrelated parallel machines

We consider the following scheduling problem. There arem parallel machines andn independent jobs. Each job is to be assigned to one of the machines. The processing of jobj on machinei requires timep _ij . The objective is to find a schedule that minimizes the makespan.Our main result is a polynomial algorithm which constructs a schedule that is guaranteed to be no longer than twice the optimum. We also present a polynomial approximation scheme for the case that the number of machines is fixed. Both approximation results are corollaries of a theorem about the relationship of a class of integer programming problems and their linear programming relaxations. In particular, we give a polynomial method to round the fractional extreme points of the linear program to integral points that nearly satisfy the constraints.In contrast to our main result, we prove that no polynomial algorithm can achieve a worst-case ratio less than 3/2 unlessP = NP. We finally obtain a complexity classification for all special cases with a fixed number of processing times.A preliminary version of this paper appeared in theProceedings of the 28th Annual IEEE Symposium on the Foundations of Computer Science (Computer Society Press of the IEEE, Washington, D.C., 1987) pp. 217–224.     

A new variable dimension simplicial algorithm for computing economic equilibria onS n ×R + m1

In this paper, a new variable-dimension simplicial algorithm is developed to compute economic equilibria on the Cartesian product of then-dimensional unit price simplexS ⁿ and them-dimensional production activity spaceR ₊ ^m . The algorithm differs from other algorithms in the number of directions in which the algorithm may leave the starting point. More precisely, the algorithm has 2 ^n+m+1 –2 rays to leave the starting point, whereas the other algorithms have at most 2 ^m (n+1) rays. The path of points generated by the algorithm can be interpreted as a globally and universally convergent price and production adjustment process. The process as well as the convergence condition are much more natural and economically meaningful than the adjustment processes obtained by other simplicial algorithms. Furthermore, we apply the algorithm to economies with linear production, economies with constant returns to scale, and economies with increasing returns to scale.This work is part of the VF-Program Competition and Cooperation.     

On parallel hashing and integer sorting

The problem of sorting n integers from a restricted range [1…m], where m is a superpolynomial in n, is considered. An o(n log n) randomized algorithm is given. Our algorithm takes O(n log log m) expected time and O(n) space. (Thus, for m = n^polylog(n) we have an O(n log log n) algorithm.) The algorithm is parallelizable. The resulting parallel algorithm achieves optimal speedup. Some features of the algorithm make us believe that it is relevant for practical applications. A result of independent interest is a parallel hashing technique. The expected construction time is logarithmic using an optimal number of processors, and searching for a value takes O(1) time in the worst case. This technique enables drastic reduction of space requirements for the price of using randomness. Applicability of the technique is demonstrated for the parallel sorting algorithm and for some parallel string matching algorithms. The parallel sorting algorithm is designed for a strong and nonstandard model of parallel computation. Efficient simulations of the strong model on a CRCW PRAM are introduced. One of the simulations even achieves optimal speedup. This is probably the first optimal speedup simulation of a certain kind.     

On the computation of multi-dimensional solution manifolds of parametrized equations

Summary A new algorithm is presented for computing vertices of a simplicial triangulation of thep-dimensional solution manifold of a parametrized equationF(x)=0, whereF is a nonlinear mapping fromR ⁿ toR ^m ,p=n–m>1. An essential part of the method is a constructive algorithm for computing moving frames on the manifold; that is, of orthonormal bases of the tangent spaces that vary smoothly with their points of contact. The triangulation algorithm uses these bases, together with a chord form of the Gauss-Newton process as corrector, to compute the desired vertices. The Jacobian matrix of the mapping is not required at all the vertices but only at the centers of certain local triangulation patches. Several numerical examples show that the method is very efficient in computing triangulations, even around singularities such as limit points and bifurcation points. This opens up new possibilities for determining the form and special features of such solution manifolds.Dedicated to Professor Ivo Babuka on the occasion of his sixtieth birthdayThis work was supported in part by the National Science Foundation under Grant DCR-8309926, the Office of Naval Research under contract N-00014-80-C-9455, and the Air Force Office of Scientific Research under Grant 84-0131     

An ODE-based approach to nonlinearly constrained minimax problems

We consider the following problem: Choosex ₁, ...,x _n to wherem ₁,m ₂,m ₃ are integers with 0m ₁m ₂m ₃, thef _i are given real numbers, and theg _i are given smooth functions. Constraints of the formg _i (x ₁, ...,x _n )=0 can also be handled without problem. Each iteration of our algorithm involves approximately solving a certain non-linear system of first-order ordinary differential equations to get a search direction for a line search and using a Newton-like approach to correct back into the feasible region when necessary. The algorithm and our Fortran implementation of it will be discussed along with some examples. Our experience to date has been that the program is more robust than any of the library routines we have tried, although it generally requires more computer time. We have found this program to be an extremely useful tool in diverse areas, including polymer rheology, computer vision, and computation of convexity-preserving rational splines.     

Parallel Implementation of Stability Analysis of Difference Schemes with MATHEMATICA

We consider a parallel algorithm for investigating the stability of the schemes of the finite-difference and finite-volume methods that approximate the two-dimensional Euler equations of compressible fluid on a curvilinear grid. The algorithm is implemented with the aid of the computer algebra system Mathematica 3.0. We apply a two-level parallelization process. At the first level, the symbolic computation of the amplification matrix is parallelized by a parallel computation of the matrix rows on different processors. At the second level, the values of the coordinates of points of the stability-region boundary are computed numerically. For the communication between the workstations, we apply a special program, LaunchSlave, which uses the MathLink communication protocol. Examples of application of the proposed parallel symbolic/numerical algorithm are presented. Bibliography: 15 titles.     

A direct method for the general solution of a system of linear equations

A computationally stable method for the general solution of a system of linear equations is given. The system isA ^Tx–B=0, where then-vectorx is unknown and then×q matrixA and theq-vectorB are known. It is assumed that the matrixA ^T and the augmented matrix [A ^T,B] are of the same rankm, wheremn, so that the system is consistent and solvable. Whenm<n, the method yields the minimum modulus solutionx _m and a symmetricn ×n matrixH _m of rankn–m, so thatx=x _m+H _my satisfies the system for ally, ann-vector. Whenm=n, the matrixH _m reduces to zero andx _m becomes the unique solution of the system.The method is also suitable for the solution of a determined system ofn linear equations. When then×n coefficient matrix is ill-conditioned, the method can produce a good solution, while the commonly used elimination method fails.This research was supported by the National Science Foundation, Grant No. GP-41158.     

A new graph triconnectivity algorithm and its parallelization

We present a new algorithm for finding the triconnected components of an undirected graph. The algorithm is based on a method of searching graphs called open ear decomposition. A parallel implementation of the algorithm on a CRCW PRAM runs inO(log² n) parallel time usingO(n+m) processors, wheren is the number of vertices andm is the number of edges in the graph.A preliminary version of this paper was presented at the19th Annual ACM Symposium on Theory of Computing, New York, NY, May 1987.Supported by NSF Grant DCR 8514961.Supported by NSF Grant ECS 8404866 and the Semiconductor Research Corporation Grant 86-12-109.     

Computing optimal preemptive schedules for parallel tasks: linear programming approaches

 We study the problem of scheduling a set of n independent parallel tasks on m processors, where in addition to the processing time there is a size associated with each task indicating that the task can be processed on any subset of processors of the given size. Based on a linear programming formulation, we propose an algorithm for computing a preemptive schedule with minimum makespan, and show that the running time of the algorithm depends polynomially on m and only linearly on n. Thus for any fixed m, an optimal preemptive schedule can be computed in O(n) time. We also present extensions of this approach to other (more general) scheduling problems with malleable tasks, due dates and maximum lateness minimization. Received: November 1999 / Accepted: November 2002 Publication online: December 19, 2002 RID="⋆" ID="⋆" This work was done while the authors were associated with the research institutes IDSIA Lugano and MPII Saarbrücken and were supported in part by the Swiss Office Fédéral de l'éducation et de la Science project n 97.0315 titled ``Platform' and by EU ESPRIT LTR Project No. 20244 (ALCOM-IT)     

Selection in Monotone Matrices and Computingkth Nearest Neighbors

Anm×nmatrix =(a_{i, j}), 1≤i≤mand 1≤j≤n, is called atotally monotonematrix if for alli₁, i₂, j₁, j₂, satisfying 1≤i₁<i₂≤m, 1≤j₁<j₂≤n.[formula]We present an[formula]time algorithm to select thekth smallest item from anm×ntotally monotone matrix for anyk≤mn. This is the first subquadratic algorithm for selecting an item from a totally monotone matrix. Our method also yields an algorithm of the same time complexity for ageneralized row-selection problemin monotone matrices. Given a setS={p₁,…, p_n} ofnpoints in convex position and a vectork={k₁,…, k_n}, we also present anO(n^4/3log^c n) algorithm to compute thek_ith nearest neighbor ofp_ifor everyi≤n; herecis an appropriate constant. This algorithm is considerably faster than the one based on a row-selection algorithm for monotone matrices. If the points ofSare arbitrary, then thek_ith nearest neighbor ofp_i, for alli≤n, can be computed in timeO(n^7/5 log^c n), which also improves upon the previously best-known result.     

Parallel Itoh–Tsujii multiplicative inversion algorithm for a special class of trinomials

In this contribution, we derive a novel parallel formulation of the standard Itoh–Tsujii algorithm for multiplicative inverse computation over the field GF(2 ^m ). The main building blocks used by our algorithm are: field multiplication, field squaring and field square root operators. It achieves its best performance when using a special class of irreducible trinomials, namely, P(x) = x ^m  + x ^k  + 1, with m and k odd numbers and when implemented in hardware platforms. Under these conditions, our experimental results show that our parallel version of the Itoh–Tsujii algorithm yields a speedup of about 30% when compared with the standard version of it. Implemented in a Virtex 3200E FPGA device, our design is able to compute multiplicative inversion over GF(2¹⁹³) after 20 clock cycles in about 0.94 μS.