mm-424 === Subject: : Re: INVERSE OF COMPLEX MATRIXUse APL2 IBM , developed at UCLA.APL has an operator (domino) who obtain directly theinverse of a matrix , real or complex.Your problem is very easy using APL.Peter Spellucci escribi.97 en el mensaje >I am looking for a routine that can find the > inverse of any matrix. Preferably in fortran . Is there aybody who >can assist me with it. > http://www.netlib.org/lapack/zgetri.f plus zgetrf.f you can use complex lu decomposition plus n solves for the unit vectors as right hand side. also the same in /linpack/zgedi.f plus zgeco.f and more hth peter=== === Subject: : VERY Fast Expression EvaluatorI've just posted an experimental math expression evaluatorwhich compiles a given expression into native machine code(Intel Pentium w/FPU). Because of the optimization strategyand management of the FPU the resulting evaluator istypically faster than the compiler generated code one wouldnormally benchmark against. Optimizations include: constantfolding, strength reduction, code inlining, no stack frames,FPU stack management, etc. A mini-disassembler is includedin the demo program so the executable code string can beexamined immediately. A DLL and demo can be found at:http://www.crbond.com/parsers.htm#MathExpr_PCE--Nobody doesn't like fast software.--Democracy: The triumph of popularity over principle.--http://www.crbond.com=== === Subject: : Plane FittingMy math has become quite rusty. I am having problems fitting a planeto a set of 3d points. To use LSF I create matrix A with its rows ofthe form [x_i-xbar y_i-ybar z_i-zbar]. Once I have found the smallesteigenvalue of matrix A using SVD, how can I find its correspondingeigenvector?C. Doe.=== === Subject: : Re: Plane FittingMy math has become quite rusty. I am having problems fitting a planeto a set of 3d points. To use LSF I create matrix A with its rows ofthe form [x_i-xbar y_i-ybar z_i-zbar]. Once I have found the smallesteigenvalue of matrix A using SVD, how can I find its correspondingeigenvector?C. Doe. let X(i) denote the ith point with coordinates x_i,y_i,z_i in your notation. then finding the minimal orthogonal least squares distance solution is equivalent to minimize the sum sum_i (A^T*X(i)-gamma)^2 subject to some normalization condition on A, with respect to the vector A and the scalar gamma, A^T*X(i) is the scalar product a_1*x_i+a2*y_i+a3*z_i hence you build a matrix with the rows X(i),1. (four columns and as much rows as you have points) then decompose by SVD this matrix A=U*S*V' where S has the same shape as A, U and V are unitary, S nonzero and nonnegative at most at the diagonal position. in your situation S should have three nonzero entries at least (the singular values of A) now you take the column of V which corresponds to the smallest singular value and normalize it such that the vector of the first three components has length one. then you got the hesse normal form of your plane A^T*x+gamma =0 length(A)=1 hth peter=== === Subject: : Re: Plane Fitting My math has become quite rusty. I am having problems fitting a plane to a set of 3d points. To use LSF I create matrix A with its rows of the form [x_i-xbar y_i-ybar z_i-zbar]. Once I have found the smallest eigenvalue of matrix A using SVD, how can I find its corresponding eigenvector?Well I'd solve (by hand probably) the eigenvalue problemAx = lambda*xwhere lambda is the smallest eigenvalue you found. Note that in theequation(A - lambda*I)x = 0the matrix A -lambda*I is singular and thus the eigenvector x can onlybe determined up to a constant.-John=== === Subject: : Re: Plane Fitting> My math has become quite rusty. I am having problems fitting a plane> to a set of 3d points. To use LSF I create matrix A with its rows of> the form [x_i-xbar y_i-ybar z_i-zbar]. Once I have found the smallest> eigenvalue of matrix A using SVD, how can I find its corresponding> eigenvector? Well I'd solve (by hand probably) the eigenvalue problem Ax = lambda*x where lambda is the smallest eigenvalue you found. Note that in the equation (A - lambda*I)x = 0 the matrix A -lambda*I is singular and thus the eigenvector x can only be determined up to a constant. -JohnYes. but I am trying to write a program to find the vector. whatmethod should I use?=== === Subject: : Re: Plane Fitting> My math has become quite rusty. I am having problems fitting a plane> to a set of 3d points. To use LSF I create matrix A with its rows of> the form [x_i-xbar y_i-ybar z_i-zbar]. Once I have found the smallest> eigenvalue of matrix A using SVD, how can I find its corresponding> eigenvector?> Well I'd solve (by hand probably) the eigenvalue problem> Ax = lambda*x> where lambda is the smallest eigenvalue you found. Note that in the> equation> (A - lambda*I)x = 0> the matrix A -lambda*I is singular and thus the eigenvector x can only> be determined up to a constant.> -John Yes. but I am trying to write a program to find the vector. what method should I use?I'd recommend PETSc from the DOD:http://www-unix.mcs.anl.gov/petsc/petsc-2/They have an extensive suite of linear solvers, as well as aneigenvalue solver I believe. I have never used it for such a purpose,but have had good success with their iterative solvers in real-worldapplications.-John=== === Subject: : Re: Plane Fitting Yes. but I am trying to write a program to find the vector. what method should I use?To some degree this depends on what the fitting targetis. If you can accept the lowest sum of squared absoluteor relative error then Gaussian least-squares by solutionof sets of linear equations is the fastest and, dependingon the solution method, most accurate and global. If youprefer fitting to the lowest sum of absolute value ofabsolute or relative error rather than SSQ, or even to thelowest peak absolute value of absolute or relative error,you'll need to use iterative non-linear solvers. Feedingthe non-linear solvers an initial parameter estimation fromthe linear SSQ fit would simple and in practice fast as well. James lips http://zunzun.com=== === Subject: : Series coeff. by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id hBUK7k424406;Hello!I want to know how we can find the coefficients of this serie? _L_ C =/_ (a_n)^2 * f_n(t) n=0C = is a constanta_n = series coefficients (we want to find)f_n = constitue an orthogonal basis functionSarah=== === Subject: : Re: Series coeff. Hello! I want to know how we can find the coefficients of this serie? _L_ C =/_ (a_n)^2 * f_n(t) n=0 C = is a constant a_n = series coefficients (we want to find) f_n = constitue an orthogonal basis functionSimply take the inner product of both sides with f_i(t) for each ifrom 0 to L. _L_ (C,f_n) =/_ (a_n)^2 * (f_n,f_i) n=0Since the {f_n} form an orthogonal basis, you have (f_n,f_i) = 0 wheni not equal to n.Thus, (C,f_n) = (a_i)^2 (f_i,f_i), from which you can determine thea_i.Dave=== === Subject: : Re: Compiling Lapack Blas with Intel Fortran (IFC) Has anyone been able to compile the BLAS and LAPACK libraries with Intel Fortan 8.0?. What were the compilation flags that you used? If this has been discussed before, could you direct me to a link/discussion thread ? -NyarlathotepYou might use Intel's Math Kernel Library (MKL) http://www.intel.com/software/products/mkl/for this purpose. It's free and optimized for Intel architectures,and it provides the same routines (I think) as the BLAS and LAPACK. If you really need to compile them by hand, can you post what errormessage you are getting?-John=== === Subject: : DeterminantsI've been fiddling round with a problem I found on mathnet for a whilenow and wondered if you guys could help:If you have a determinant | 1 1 1 ||Sin (x+a) Sin (x+b) Sin (x+c) ||Cos (x+a) Cos (x+b) Sin (x+c) |How can you show it is independent of x?SW=== === Subject: : Re: Determinants I've been fiddling round with a problem I found on mathnet for a while now and wondered if you guys could help: If you have a determinant | 1 1 1 | |Sin (x+a) Sin (x+b) Sin (x+c) | |Cos (x+a) Cos (x+b) Sin (x+c) | How can you show it is independent of x? SWAs you have written it, the expanded determinant is sinC ( sin B - cosB - sinA + cosA) - sinA cosB + cosA sinBwhere A = x+a, B = x+b C = x+c . The last 2 terms are -sin(a-b)which is independent of x. However, as you can see by differentiatingw.r.t. x, the first 4 terms are not independent of x.However, if, as the first respondent said, you should have written | 1 1 1 | |Sin (x+a) Sin (x+b) Sin (x+c) | |Cos (x+a) Cos (x+b) Cos (x+c) | ^^^^^^^^^the result is independent of x. ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/=== === Subject: : Re: Determinants I've been fiddling round with a problem I found on mathnet for a while now and wondered if you guys could help: If you have a determinant | 1 1 1 | |Sin (x+a) Sin (x+b) Sin (x+c) | |Cos (x+a) Cos (x+b) Sin (x+c) |That last entry should be cos(x+c) not sin(x+c), right? How can you show it is independent of x?Hint: The determinant of a matrix of the form (1,1,1; a,b,c; d,e,f)is related to the area of the triangle with vertices (a,d) etc.Alternative hint (don't try to use them both at once): brute forceusing trigonometric identities.=== === Subject: : Re: Determinants> I've been fiddling round with a problem I found on mathnet for a while> now and wondered if you guys could help:> If you have a determinant > | 1 1 1 |> |Sin (x+a) Sin (x+b) Sin (x+c) |> |Cos (x+a) Cos (x+b) Sin (x+c) | That last entry should be cos(x+c) not sin(x+c), right?> How can you show it is independent of x? Hint: The determinant of a matrix of the form (1,1,1; a,b,c; d,e,f) is related to the area of the triangle with vertices (a,d) etc. Alternative hint (don't try to use them both at once): brute force using trigonometric identities.I think this operation may help ... = ( sin(x+b)cos(x+c) - sin(x+c)*cos(x+b) ) - ( sin(x+a)cos(x+c) -sin(x+a)*cos(x+c) ) +(sin(x+a)cos(x+b) - sin(x+a)*cos(x+b)) = sin((x+b) - (x+c)) - sin((x+a) - (x+c)) + sin((x+a) - (x+b)) =sin(b-c)-sin(a-c)+sin(a-b)we used sin(z-y) = sin(z)cos(y) - sin(y)cos(z).=== === Subject: : Re: Determinantsyes, it is sin(x+c).I've expanded the determinant fully and got a big mess with none of myterms cancelling or using sin^2+cos^2=1. I don't see how your triangleidea helps. I know the determinant is twice that of the correspondingtriangle but it doesn't show any independence of x.Another hint would be very helpful!SW> I've been fiddling round with a problem I found on mathnet for a while> now and wondered if you guys could help:> If you have a determinant > | 1 1 1 |> |Sin (x+a) Sin (x+b) Sin (x+c) |> |Cos (x+a) Cos (x+b) Sin (x+c) | That last entry should be cos(x+c) not sin(x+c), right?> How can you show it is independent of x? Hint: The determinant of a matrix of the form (1,1,1; a,b,c; d,e,f) is related to the area of the triangle with vertices (a,d) etc. Alternative hint (don't try to use them both at once): brute force using trigonometric identities.=== === Subject: : Re: Determinantsyes, it is sin(x+c).I've expanded the determinant fully and got a big mess with none of myterms cancelling or using sin^2+cos^2=1. I don't see how your triangleidea helps. I know the determinant is twice that of the correspondingtriangle but it doesn't show any independence of x.Another hint would be very helpful!SW>I've been fiddling round with a problem I found on mathnet for a while>now and wondered if you guys could help:>If you have a determinant >| 1 1 1 |>|Sin (x+a) Sin (x+b) Sin (x+c) |>|Cos (x+a) Cos (x+b) Sin (x+c) |> That last entry should be cos(x+c) not sin(x+c), right?>How can you show it is independent of x?> Hint: The determinant of a matrix of the form (1,1,1; a,b,c; d,e,f)> is related to the area of the triangle with vertices (a,d) etc.> Alternative hint (don't try to use them both at once): brute force> using trigonometric identities. Expand by minors of the first row. If Air France must cancel Paris-LAX Christmas Eve flights, let everyone ride with Santa Claus. He can stop at his North Pole home en route. pmontgom@cwi.nl Home: San Rafael, California Microsoft Research and CWI=== === Subject: : Re: Determinants yes, it is sin(x+c).Um. Do you mean yes, it is cos(x+c) or no, it is sin(x+c)?!(I claimed it ought to be cos(x+c); you'd written sin(x+c).) I've expanded the determinant fully and got a big mess with none of my terms cancelling or using sin^2+cos^2=1.The big mess should include some bits on which you canuse the trigonometric sum and difference formulas.Use them backwards, I mean: spot things that looklike sin(x)cos(y)-sin(y)cos(x) or whatever. I don't see how your triangle idea helps. I know the determinant is twice that of the corresponding triangle but it doesn't show any independence of x.Oh yes it does :-). You have three points like (cos(A),sin(A)) forthree different values of A, right? So they form a triangleinscribed in the unit circle. What happens to that trianglewhen you change the value of x?=== === Subject: : Re: Determinants> yes, it is sin(x+c). Um. Do you mean yes, it is cos(x+c) or no, it is sin(x+c)?! (I claimed it ought to be cos(x+c); you'd written sin(x+c).)> I've expanded the determinant fully and got a big mess with none of my> terms cancelling or using sin^2+cos^2=1. The big mess should include some bits on which you can use the trigonometric sum and difference formulas. Use them backwards, I mean: spot things that look like sin(x)cos(y)-sin(y)cos(x) or whatever.Another idea is to forget about the trigonometric at this stage.Instead, consider the det to be some function of x, i.e. L(x) = det [M(x)]Now instead of trying to prove that L(x) = constant, try taking derivatives.It is almost trivial to show that dL/dx = 0, which is just another way ofsayingthat det[M(x)] is independent of x.Personally, I prefer Gareth's triangle method to an algebraic route.~Glynne> I don't see how your triangle> idea helps. I know the determinant is twice that of the corresponding> triangle but it doesn't show any independence of x. Oh yes it does :-). You have three points like (cos(A),sin(A)) for three different values of A, right? So they form a triangle inscribed in the unit circle. What happens to that triangle when you change the value of x? -- Gareth McCaughan .sig under construc=== === Subject: : Partial differentiation and derivatives.Can any of you help me with this problem I've been having?If you have a function v(r)=e^(k(r-x)) in which r^2=x^2+y^2 is thereany way you can find the partial derivative of dv/dy at constant xwithout substitution. I found dv/dx as -ke^k((r-x)) but did it withthe chain rule.(^ is to the power of)Jon=== === Subject: : Re: Partial differentiation and derivatives. Can any of you help me with this problem I've been having? If you have a function v(r)=e^(k(r-x)) in which r^2=x^2+y^2 is there any way you can find the partial derivative of dv/dy at constant x without substitution. I found dv/dx as -ke^k((r-x)) but did it with the chain rule. Sadly, you did it wrong. That's not the answer. v depends on x in TWO places. The chanin rule is the right thing to use. dv/dx = k* (x/r - 1) * v(r,x) Hence dv/dy = k * (y/r) * v(r,x) ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/=== === Subject: : 2d interpolationDear All, Could anyone give me the reference about 2d interpolation? It is betterthat it contains the fast algorithm. ThxTim=== === Subject: : Re: 2d interpolation Could anyone give me the reference about 2d interpolation? It is better that it contains the fast algorithm.Look at the data fitting section in http://www.mat.univie.ac.at/~neum/stat.html[~ is a tilde]Arnold Neumaier=== === Subject: : Re: 2d interpolation >Dear All, > Could anyone give me the reference about 2d interpolation? It is better >that it contains the fast algorithm. Thx >Tim >what kind of interpolation? what kind of data? rectangular grid? scattered data? smooth piecewise polynomial?some proposals:In Fortran 77, which should be f90 compatible:http://www.netlib.org/toms/684C1 and C2 Interpolation on Triangles with Quintics and Nonics.there is more in netlib/toms , e.g. piecewise linear and continuous.on rectangular grids, you may have easier solutions.hthpeter=== === Subject: : Re: 2d interpolationNumerical Recipes (cf NAG: Cambridge U.)~ Includes C::Fortan::PascalBryan Dear All, Could anyone give me the reference about 2d interpolation? It is better that it contains the fast algorithm. Thx Tim=== === Subject: : Re: 2d interpolationThx Numerical Recipes (cf NAG: Cambridge U.)~ Includes C::Fortan::Pascal Bryan> Dear All,> Could anyone give me the reference about 2d interpolation? It isbetter> that it contains the fast algorithm. Thx> Tim=== === Subject: : Conjugate Gradients on Poisson's equation - spectral error analysisI have a Poisson equation (Ax=b) that arises from fluid dynamics. I'musing the conjugate gradient (CG) method to solve for x.Now, let's discuss two alternative techniques:1.) Gauss-(Seidel) (or SO) relaxation. This technique is nice forsolving high-frequency components on the grid, but sucks at at lowerfrequency components, because it takes alot of time for theinformation to propagate through the grid by the massive number ofrelaxation iterations.2.) Multigrids. As an opposite of relaxation, multigrids has gooderror properties on lower-frequency components, but does more badly onhigher-frequency components (and as a consequence, relaxation issometimes used to refine the higher-frequency components in amultigrid solution). I haven't implemented multigrids myself, but thisis what I've heard of :)So, we know the pros and cons of relaxation and multigrids for theerror at a spectral level. But what about CG, then? Does CG do betterat higher, or at lower frequency components? Or does it do equallywell on both? Should I apply a bit of relaxation after the CGiterations to ensure high-frequency accuracy?Any feedback appreciated, - Mikko Kauppila=== === Subject: : Re: Conjugate Gradients on Poisson's equation - spectral error analysis >I have a Poisson equation (Ax=b) that arises from fluid dynamics. I'm >using the conjugate gradient (CG) method to solve for x. >Now, let's discuss two alternative techniques: >1.) Gauss-(Seidel) (or SO) relaxation. This technique is nice for >solving high-frequency components on the grid, but sucks at at lower >frequency components, because it takes alot of time for the >information to propagate through the grid by the massive number of >relaxation iterations. > >2.) Multigrids. As an opposite of relaxation, multigrids has good >error properties on lower-frequency components, but does more badly on >higher-frequency components (and as a consequence, relaxation is >sometimes used to refine the higher-frequency components in a >multigrid solution). I haven't implemented multigrids myself, but this >is what I've heard of :)by the correct use of the multigrid feature it should smooth components at allfrequencies fast, the problem being the construction of the nested grid? >So, we know the pros and cons of relaxation and multigrids for the >error at a spectral level. But what about CG, then? Does CG do better >at higher, or at lower frequency components? Or does it do equally >well on both? Should I apply a bit of relaxation after the CG >iterations to ensure high-frequency accuracy? > >Any feedback appreciated, > - Mikko Kauppilathe error behaviour of CG depends strongly on the eigenvalue distribution of the matrix and for the error in x_k there holds (x_k-x*)^TA(x_k-x*) = min (x_0-x*)^TA(I+AP_{k-1}(A))^2(x_0-x*) <= min max_i |1+lambda_iP_{k-1}(lambda_i))^2 (x_0-x*)^TA(x_0-x*) the minimum being taken over all polynomials of degree k-1.and lambda_i the i-th eigenvalue of A Lambda_1>=...>=lambda_n>0.unfortunately the matrix for the discrete poisson problem has no goodeigenvalue distribution in this sense and for this reason you should apply cgwith preconditioning only. with a good preconditioner cg will be competitivewith multigridhthpeter=== === Subject: : Re: Conjugate Gradients on Poisson's equation - spectral error analysis 1.) Gauss-(Seidel) (or SO) relaxation. This technique is nice for solving high-frequency components on the grid, but sucks at at lower frequency components, because it takes alot of time for the information to propagate through the grid by the massive number of relaxation iterations.Sounds correct. 2.) Multigrids. As an opposite of relaxation, multigrids has good error properties on lower-frequency components, but does more badly on higher-frequency components (and as a consequence, relaxation is sometimes used to refine the higher-frequency components in a multigrid solution). I haven't implemented multigrids myself, but this is what I've heard of :)Relaxation is *always* used. So I guess MG will be pretty good at highfrequencies too. So, we know the pros and cons of relaxation and multigrids for the error at a spectral level. But what about CG, then? Does CG do better at higher, or at lower frequency components? Not sure about the full answer, but CG builds up an approximation of thespectrum that converges quickest in the outer eigenvalues. I guess thatmeans the middle is the problem. Should I apply a bit of relaxation after the CG iterations to ensure high-frequency accuracy?Won't hurt. Not sure how much good it will do.V.email: lastname at cs utk eduhomepage: cs utk edu tilde lastname=== === Subject: : Round off errorDistribution: inetHi Folks:I was wondering if there is some sort of test that can be run on anumerical algorithm to quantify round off error. Needless to say, thealgorithm is stable.For me, this comes up as follows: I run finer and finer meshes on aproblem and observe some particular value of the solution. If the exactvalue is y_ex, the computed values approach and then diverge from y_ex.I am wondering if I can attribute this divergence to round off errorsentirely. Is there some other definitive test that will establish this ?your help.Suresh--It's all there when you look,Into Santa's book.=== === Subject: : Re: Round off error Hi Folks: I was wondering if there is some sort of test that can be run on a numerical algorithm to quantify round off error. Needless to say, the algorithm is stable. For me, this comes up as follows: I run finer and finer meshes on a problem and observe some particular value of the solution. If the exact value is y_ex, the computed values approach and then diverge from y_ex. I am wondering if I can attribute this divergence to round off errors entirely. Is there some other definitive test that will establish this ? your help. SureshWhy don't you use Interval Arithmetic. In that case you will alwayshave an upper bound on the round off errors in your comnputation.Henrik=== === Subject: : Re: Round off errorDistribution: inet >Hi Folks: >I was wondering if there is some sort of test that can be run on a >numerical algorithm to quantify round off error. Needless to say, the >algorithm is stable. >For me, this comes up as follows: I run finer and finer meshes on a >problem and observe some particular value of the solution. If the exact >value is y_ex, the computed values approach and then diverge from y_ex. >I am wondering if I can attribute this divergence to round off errors >entirely. Is there some other definitive test that will establish this ? >your help. >yes, with no universal answer. if you use interval arithmetic instead of normal rounded one, you get true inclusions of the error, but often the bounds become too pesimistic. if you know testcases withexact answers, you can get a first impression by comparing the resluts withthe exact ones (using exact inputs of course). but the result will be often to optimistic. one brute force approach is to use the same algorithm on the samedata using different precisions, again on a (carefully) selected set of testcases. what you observe is the normal effect of accumulation of roundoff errors. for example, if you have to compute a derivative by finite differences, then depending on the order of the formula you use you will obeservethis effect (for the ordinary simple forward difference an optimal stepsize (up to a constant) is sqrt(eps) with eps describing the precision of thecomputation. if you have to solve an initial value problem of an ODE you have the same effect: roundoff influence grows as N*eps, N the step number, whereas discretization error goes down as 1/N^p, p the consistency order, and again there is an optimal N (an optimal stepsize).hthpeter=== === Subject: : Re: Round off errorHi Folks:I was wondering if there is some sort of test that can be run on anumerical algorithm to quantify round off error. Needless to say, thealgorithm is stable.For me, this comes up as follows: I run finer and finer meshes on aproblem and observe some particular value of the solution. If the exactvalue is y_ex, the computed values approach and then diverge from y_ex.I am wondering if I can attribute this divergence to round off errorsentirely. Is there some other definitive test that will establish this ?your help.SureshCan you run the code in a higher precision - it is sort of cheatingbut it is one way to check these things out. If you are runningdouble precision on an intel machine you could try going to extendedprecision. You might even play with reducing the precision to singleprecision to see if that shows how sensitive the calculation is to theprecision of the calculation.Is there a proof that the method you are using is stable and willconverge? What is the order of convergence? Can you compare theerrors as you change the mesh with the expected rate of convergence.Basically it seems to me that if the numerical technique has by arigorous proof been shown to be stable and to converge then thebehaviour you are describing must be the result of either:1. round off error2. the problem you are trying to solve is ill conditioned orsingular.=== === Subject: : Re: Round off error Hi Folks: I was wondering if there is some sort of test that can be run on a numerical algorithm to quantify round off error. Needless to say, the algorithm is stable. For me, this comes up as follows: I run finer and finer meshes on a problem and observe some particular value of the solution. If the exact value is y_ex, the computed values approach and then diverge from y_ex. I am wondering if I can attribute this divergence to round off errors entirely. Is there some other definitive test that will establish this ? your help. Suresh -- It's all there when you look, Into Santa's book.Are you, by any chance, fitting data with least-squares? If so, you canexpect that behavior. ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ God is not willing to do everything, and thereby take away our free wiil and that share of glory that rightfully belongs to ourselves. -- N. Machiavelli, The Prince.=== === Subject: : Re: another sparse eigenvector question I'm writing a paper, and I don't want to say anything wrong. Given a sparse symmetric matrix M (M is n by n and has O(n) non-zero entries), in the worst case is it possible that a Cholesky decomposition of M will require O(n^2) memory and O(n^3) time? Yes: take a matrix with full first row, first column, and diagonal, and otherwise empty. I understand that it will generally be pretty good in practice if you reorder the matrix appropriately.Yes, but things depend a lot on the structure of the sparsity graph.Arnold Neumaier=== === Subject: : virtual displacmentsWhat is a virtual displacement as appeased to a functional displacement=== === Subject: : statistics/regression experimentHiya,I'm doing some statistics and have got into a rut. I'd reallyappreciate any help that any of you can give me.If you have a set of points from an experiment e.g. (x1,y1), (x2,y2),(x3,y3)....(xn,yn) how do you find the best straight line through thepoints? I.e. how do you find p and q so that the sum of(yi-p*(xi)-q)^2 is a minimum. The sum if from i=1 to n. (i is supposedto be in subscript).Jon=== === Subject: : Re: statistics/regression experimentI'm doing some statistics and have got into a rut. I'd reallyappreciate any help that any of you can give me.If you have a set of points from an experiment e.g. (x1,y1), (x2,y2),(x3,y3)....(xn,yn) how do you find the best straight line through thepoints? I.e. how do you find p and q so that the sum of(yi-p*(xi)-q)^2 is a minimum. The sum if from i=1 to n. (i is supposedto be in subscript).Jon===========================log on to:http://people.hofstra.edu/faculty/Stefan_Waner/newgraph/ regressionframes.htmlDick=== === Subject: : Re: how to design iterative optimization algorithms to solve matrix equation? Any books?>Dear all,>I am facing with this difficult problem. Please help me!>In this problem, I need to construct some matrices which satisfy the>following matrix equation:>(( A * A )V1 + ( B * B) V2) V = (D * D)>where * denotes Kronecker product, A, B, V1, V2, V are unknownmatrices>that needs solving; they are all square. Some structure needs to beimposed:>I have certain pattern for A and B; and V1, V2, and V are required to be>diagonal... Matrix D is given...>The task is to find the best approximation of the above-mentioned A, B,V1,>V2 and V... I have been thinking about this for long time... can anybody>give me some hints?>I guess it is difficult to find closed form analytical solution.. How diificult is it - it seems to me that the problem has an awful lot of structure.>. How to>design iterative algorithm to let computer search for the answer? It is>really hard... please help me! You are trying to solve the nonlinear system of equations (( A * A )V1 + ( B * B) V2) V - (D * D) = 0 If the problem is not too large use Newton's method. If the problem is too big for that look for a canned program that will solve nonlinear systems without needing to generate and then solve the jacobian.Dear Joe,Thank you for your answer. However I still got questions:a) (( A * A )V1 + ( B * B) V2) V - (D * D) = 0These * are Kronecker products not the normal multiplier, how to designNewton's iteration algorithm for this equation?b) The worst thing is that this problem is not continuous; the elements in Vcan be arbitary; but the elements in A, B, V1, V2 need to be integer...hence I lost completely how to do iterative algorithm for these kind ofmix-integer-continuous problem...Could you please help me more?-Walalla=== === Subject: : what is typical good algorithm to attack integer minimization problems?Dear all,My problem is to construct some integer matrices of certain pattern tosatisfy a matrix equation... There are so many unknowns, and there are alsoso many equations, (number of equations > number of unknowns)...Since the required matrices are integer, I cannot use Newton, or Conjugatedescent methods, etc... I have to do some search algorithm... but sincethere are so many unknowns, the search space is too large...So I am thinking of having an initial matrix, then do some change, see ifthe matrix equation(left hand side - right hand side) approximately morenear zero and find the minimization point...But even this, still make a large search space when the size of matrix getslarger... moreover, another constraint is that the number of non-zeroelements of this matrix should be as small as possible...Please tell me what kind of problem is this? I even don't know how toPlease give me some pointers!-Walala=== === Subject: : Re: what is typical good algorithm to attack integer minimization problems? Dear all, My problem is to construct some integer matrices of certain pattern to satisfy a matrix equation... There are so many unknowns, and there are also so many equations, (number of equations > number of unknowns)... Since the required matrices are integer, I cannot use Newton, or Conjugate descent methods, etc... I have to do some search algorithm... but since there are so many unknowns, the search space is too large... So I am thinking of having an initial matrix, then do some change, see if the matrix equation(left hand side - right hand side) approximately more near zero and find the minimization point... But even this, still make a large search space when the size of matrix gets larger... moreover, another constraint is that the number of non-zero elements of this matrix should be as small as possible... Please tell me what kind of problem is this? I even don't know how toPossible keywords:(linear or lonlinear) mixed integer programmingconstraint programmingNEOSWriting your own algorithm is likely to be inefficient unless you are very experienced. Thus take one of the packages that already exist!Arnold Neumaier=== === Subject: : Re: what is typical good algorithm to attack integer minimization problems?Originator: nmm1@virgo.cus.cam.ac.uk||> My problem is to construct some integer matrices of certain pattern to|> satisfy a matrix equation... There are so many unknowns, and there are also|> so many equations, (number of equations > number of unknowns)...||> Since the required matrices are integer, I cannot use Newton, or Conjugate|> descent methods, etc... I have to do some search algorithm... but since|> there are so many unknowns, the search space is too large...||> So I am thinking of having an initial matrix, then do some change, see if|> the matrix equation(left hand side - right hand side) approximately more|> near zero and find the minimization point...||> But even this, still make a large search space when the size of matrix gets|> larger... moreover, another constraint is that the number of non-zero|> elements of this matrix should be as small as possible...||> Please tell me what kind of problem is this? I even don't know how to||> Possible keywords:||> (linear or lonlinear) mixed integer programming||> constraint programmingSee Knuth, volume II, for the Spectral test. It does preciselythis. There are a zillion other potential keywords, too, becauseit is a common problem. For example, much of seriation and evenmulti-dimensional scaling uses these techniques and has done sincethe 1950s. Another relevant concept is linear programming.|> NEOS||> Writing your own algorithm is likely to be inefficient unless |> you are very experienced. Thus take one of the packages that |> already exist!Grrk. Well, yes, but ....This is a classically foul problem, because generic solutions justdon't work. As a result, virtually EVERY useful approach is a hack,and the difficulty is getting one that will work on YOUR problem.You often need to roll your own, whether or not you have the skillsto do it - that is, after all, how most people get into this area!But I agree that trying the standard ones BEFORE starting to hackis a good idea, whether you are experienced or not. I may be a madhacker from way back, but even I will take a look around for existingcode first :-)Nick Maclaren.=== === Subject: : Re: Volume Analysis?This method below forces the ellipsoid to be aligned with the principal axesand won't necessarily provide the wished *enclosed volume*!Perhaps a rotation of the data points along their inertia principal axesbefore applying the method would provide a reasonably optimal ellipsoidof the average data point cloud, but this is not the wished result.surface needs to be convex or not, and differentiable or not.If not necessarily convex and differentiable Voronoi decomposition providesan optimal result. If the surface may be a polyhedron but must beIf it must be the smallest enclosing sphere or ellipsoid then a non-linearminimisation technique may be applied. Dan>I have a set of 3D data points. I want to fit a surface to these>points and calculate the enclosed volume. For example, find the best>fit ellipsoid or sphere.>Can someone point me in the direction of what algorithims are used and>better yet, if matlab has any functions that can do this?>Please just post here so everyone can benefit from the response. I would probably use least squares fitting. Write the position of a point on an ellipsoid as x = x_0 + a sin ( theta ) cos ( phi ) y = y_0 + b sin ( theta ) sin ( phi ) z = z_0 + c cos ( theta ) Now for each experimental 3-dim point vec r _k calculate theta _k and phi _k , and then compute (and minimize) chi ^2 = sumlimits_{k = 1}^N {left( {frac{{left[ {vec rleft( {theta _k ,phi _k } right) - vec r_k } right]^2 }}{{sigma _k^2 }}} right)} where the sigma _k^2 is the estimated variance of each experimental point vec r_k . You have 6 free parameters that you are varying in the minimization routine: the position of the geocenter, x_0, y_0, z_0, and the 3 radial parameters, a, b, c. Since the fitting function vec r ( theta , phi ) is linear in these parameters, the problem is just linear least-squares. Once you have minimized, you can solve the 6 linear equations for the unknown parameters. You only need a,b,c since the volume of such an ellipsoid is just V = frac{{4 pi} {3}} abc .=== === Subject: : Re: Approximate GCD, any package? I am looking for software packages for computing approximate GCD of polynomials whose coefficients may be inexact. I'd appreciate any Multivariate approximate GCD-finders are of particular interest. I found EpsilonGCD, QuasiGCD and QRGCD in Maple SNAP package for univariate polynomials only. Scilab function gcd can handle inexact data, also for the univariate case. Exact GCD for multivariate polynomials are quite standard in computer algebra packages. What if the coefficient are approximately known?The problem becomes ill-posed, unless you have bounds on the coefficients.If you know the expected degree of the gcd, you can phrase the factorization as an overdetermined nonlinear system of equations for the coefficients and try to solve it as a nonlinear least squares problem. If you have bounds on the coefficients, you cantry to solve it instead as a constrained nonlinear program with trivial objective function. In both cases you may miss the solution due to spurious local minima - in that case you'd need a global solver such as BARON.If you don't know the degree you need to do this for all resonable degrees.You might also look at papers by Stetter who analyzed exactcomputer algebra algorithms in the presence of uncertainties- his work might contain better algorithms than what I outlined above.Arnold Neumaier=== === Subject: : Re: Approximate GCD, any package?The more proper way to say may be: GCD-finding is ill-posed if it is not well-posed.Knowing the bounds on the coefficients or not may not change that.You are absolutely right that GCD-finding is essentially a nonlinear leastsquares problem. However, you don't NEED to know the degree of GCD. All theGCD information is in the coefficients, you just need to figure it out by aproper algorithm.A similar case is polynomial root-finding. Many believed that you need priorknowledge on the multiplicity before you calculate the roots right. In fact,the multiplicity information is hiding in the coefficients and can be diggedout computationally.I actually have two methods and testing codes for computing approximate GCD.One for univariate polynomials and the other is for multivariate cases. Thereason I posted the question is to see if I miss someone else's software tocompare with.Zeke The problem becomes ill-posed, unless you have bounds on the coefficients. If you know the expected degree of the gcd, you can phrase the factorization as an overdetermined nonlinear system of equations for the coefficients and try to solve it as a nonlinear least squares problem. If you have bounds on the coefficients,you can try to solve it instead as a constrained nonlinear program with trivial objective function. In both cases you may miss the solution due to spurious local minima - in that case you'd need a global solver such as BARON. If you don't know the degree you need to do this for all resonable degrees. You might also look at papers by Stetter who analyzed exact computer algebra algorithms in the presence of uncertainties - his work might contain better algorithms than what I outlined above. Arnold Neumaier=== === Subject: : Re: Matrix inequality/interpolation result I have asked Rob Stevenson from Utrecht and he refered me to Lemma 6.2.5 of Hackbusch's multigrid book where one can find a more general interpolation result. Lemma 6.2.5: Let $Lambda_1, Lambda_2$ be spd and $A$ be an arbitrary matrix. Then $$ norm{Lambda_2^beta A Lambda_1^{-beta}} leq norm{Lambda_2^alpha A Lambda_1^{-alpha}}^{frac{gamma-beta}{gamma-alpha}} norm{Lambda_2^gamma A Lambda_1^{-gamma}}^{frac{beta-alpha}{gamma-alpha}} $$ for $alpha leq beta leq gamma$. Proof: $log g$ is convex for $g(s)=norm{Lambda_2^s A Lambda_1^{-s}}$. This follows from $$ g(s)^2 = norm{Lambda_2^s A Lambda_1^{-2s} A^* Lambda_2^s} = rho(Lambda_2^s A Lambda_1^{-2s} A^* Lambda_2^s) leq g(r) g(t) $$ for $s=(r+t)/2$. So there is a very nice and easy proof of my question (choose $A=Id$ in this lemma).But this reduces to your problem only when the lambda's commute!?Arnold Neumaier=== === Subject: : Re: Simple sum question If I state, sum^{N}_{i=1} t_{i} = 0 what does the following sum equal? sum^{N}_{i=1} t^{2}_{i} My naive thought was that it was also zero but after working backwards from the next equation in a stats/astronomy paper, the author has something like: sum^{N}_{i=1} t^{2}_{i} = (Delta t^2)/12 N^3 Any ideas?Here's an idea: it is an approximation, after it should be O(N^2), andDelta t:=max|t_i - t_{i-1}|. This is true at least for the case of linearlyspaced t_i 's (proof below), so I guess this approx is valid if the t_i areapproximately linearly spaced.I think otherwise the author's result does not make sense. Unless, ofcourse, if you _define_ the Delta t by the required result, and that wouldmake it rather trivial.The proof in linear case: suppose t_i are linearly spaced, say t_i=a*(i-M)where M=N/2 or smth like that to make sum t_i = 0, then Delta t = a andsum t_i^2 = Delta t^2 sum (i-M)^2.Now sum (i-M)^2 = sum i^2 - 2*M*sum i + sum M^2.sum i^2 = (1/3)*N^3 +O(N^2),sum i = (1/2)*N^2 +O(N),2*M = N +O(1),sum M^2 = (1/4)*N^3 + O(N^2)now collect the constants, you get (1/12)*N^3 + O(N^2) as claimed.Cheers,Teijo.=== === Subject: : Re: Simple sum question Hello all, I am wishing that I paid more attention in my pre-calculus/calculus courses when it came to sums and series... If I state, sum^{N}_{i=1} t_{i} = 0 what does the following sum equal? sum^{N}_{i=1} t^{2}_{i} My naive thought was that it was also zero but after working backwards from the next equation in a stats/astronomy paper, the author has something like: sum^{N}_{i=1} t^{2}_{i} = (Delta t^2)/12 N^3 Any ideas?What is Delta in this case?If we let t_i = (-1)^i, and N is even we have:sum^{N}_{i=1} t_i = -1 + 1 -1 + 1 ... = 0however,sum^{N}_{i=1} t_i^2 = 1 + 1 + 1 + 1 = Nwhich provides a contradiction to the first naive thought. I haveworked with it a little but, but it is difficult to get a simplerrepresentation than t_i^2. I'm not sure how to obtain the author'sresult in any case.-John=== === Subject: : Re: Simple sum question Hello all, I am wishing that I paid more attention in my pre-calculus/calculus courses when it came to sums and series... If I state, sum^{N}_{i=1} t_{i} = 0 what does the following sum equal? sum^{N}_{i=1} t^{2}_{i} My naive thought was that it was also zero but after working backwards from the next equation in a stats/astronomy paper, the author has something like: sum^{N}_{i=1} t^{2}_{i} = (Delta t^2)/12 N^3 Any ideas?Of course. Look at the sum (LaTeX) of alternating 1's and -1's:[ sumlimits_{k = 0}^N {left( { - 1} right)^k } = left{ {begin{array}{*{20}c} {0,quad N,{rm odd}} hfill {1,quad N,{rm even}} hfill end{array}} right.]Now look at the sum of the squares of the terms (all 1's): the answer is N+1. ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ God is not willing to do everything, and thereby take away our free wiil and that share of glory that rightfully belongs to ourselves. -- N. Machiavelli, The Prince.=== === Subject: : DifferentiationGiven that f(x) <= c, then df(x)/dx <= 0. Generally, what does this mean?=== === Subject: : Re: DifferentiationIn sci.math.num-analysis, Wumin Given that f(x) <= c, then df(x)/dx <= 0. Generally, what does this mean?I'm not sure what it means, as the equationf(x) = c - x^2clearly has f(x) <= c for every x in the reals and yetf'(x) = df(x)/x = -2xis greater than 0 for negative x, less than 0 for positive x.#191, ewill3@earthlink.netIt's still legal to go .sigless.=== === Subject: : Re: Differentiation Given that f(x) <= c, then df(x)/dx <= 0. Generally, what does this mean?It isn't true. That is, the second part doesn't follow from the first. ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/=== === Subject: : Re: Differentiation Given that f(x) <= c, then df(x)/dx <= 0. Generally, what does this mean?Which part of it?The question is a strange-sounding one, and it seems possibleto me that there's something earlier that you're not quiteunderstanding right. It might be helpful to go back a step ortwo. Where does your question come from? What is it that you'retrying to understand or solve?=== === Subject: : VarianceI have a set of sample data which I wanted to compute variance ofactivity intensity over all events. But the problem is that each ofthese events was bursty in nature and has long interval of inactivity.Most events were short and disappeared before the sampling wascompleted. Furthermore, most events started at different time points.I did the following. Since the events were bursty, I chopped theevents into fixed size segments along the time horizon. The segmentsize was chosen such that it was a constant multiply of the smallestburst interval observed. I then computed the mean activity intensityof each segment for every event. The reason for such segment size wasbecause I did not wish the long inactivity interval to affect theaverage intensity. For each segment I partitioned it into smallerchunks and computed their respective average intensity.The pair-wise sub-covariance is subcov(A,B) = 1/M * sum_i {(A_i - mean(A)) * (B_i - mean(B))}where `i' is index of chunk and `M' is total number of chunks in asegment. The covariance between pair of events isCOV(X,Y) = 1/N * sum_j {subcov(X^j,Y^j)}where `j' is index of segment and N is total number of segments. Letthe set of sample data be Z. ThenVAR(Z) = sum_u sum_v COV(Z^u,Z^v)where `u' and `v' are indices of event. Note that some events startedlater and some stopped early. In this case, data in their respectivesegments where they were absence were set to zeros.I am quite unsure whether the procedure I've taken was correct. Iwould like to ask for advice.Thank you in advance.=== === Subject: : DiffGiven that f(x) <= c, then df(x)/dx <= 0. Generally, what does this mean?=== === Subject: : Re: Diff Given that f(x) <= c, then df(x)/dx <= 0. Generally, what does this mean?It means you should do your own homework.=== === Subject: : SortingAny pointers to newsgroups/forums that entertain serious questions aboutsorting?Heck, Numerical Recopies in C mentions sorting, so I thought that I'd askhere.http://www.standards.com/; See Howard Kaikow's web site.=== === Subject: : Re: Sorting Any pointers to newsgroups/forums that entertain serious questions about sorting? Heck, Numerical Recopies in C mentions sorting, so I thought that I'd ask here. -- http://www.standards.com/; See Howard Kaikow's web site.Why not ask your questions _here_ ?And (as always) read:Donald E. Knuth: The Art of Computer Programming, Volume 3 /Sorting and Searching. Second Edition. Addison-WesleyHugo Pfoertner=== === Subject: : Re: SortingOne of my quests is to find a source, not code, that explains the differencebetween a ripple sort and a bubble sort.I have Sedgewick's book, no mention of ripple sort in the index.No mention of ripple sort in the index of the Knuth edition I have.I have other, more interesting quests, but wish to resolve this one first.http://www.standards.com/; See Howard Kaikow's web site.> Any pointers to newsgroups/forums that entertain serious questions about> sorting?> Heck, Numerical Recopies in C mentions sorting, so I thought that I'dask> here.> --> http://www.standards.com/; See Howard Kaikow's web site. Why not ask your questions _here_ ? And (as always) read: Donald E. Knuth: The Art of Computer Programming, Volume 3 / Sorting and Searching. Second Edition. Addison-Wesley Hugo Pfoertner=== === Subject: : Re: Sorting<I dont know any book, but the ripple sort takes the first element compares it with each of the remaining and swap them when the one in higher position is greater. At the end of the process the first element will be the lowest, so the algoritm repeats starting with the second.The bubble sort compares the elements pairwise (fist w/ second then second w/ thisr) swap them if the one in higher position is greater. After the first step the highest element is at the end. The cycle is repeated until there is no change during and cycle.=== === Subject: : Re: Sorting One of my quests is to find a source, not code, that explains the difference between a ripple sort and a bubble sort. I have Sedgewick's book, no mention of ripple sort in the index. No mention of ripple sort in the index of the Knuth edition I have. I have other, more interesting quests, but wish to resolve this one first. -- http://www.standards.com/; See Howard Kaikow's web site.[...]You can find a lot of material in Google Groups, e.g. a discussionAFAIK ripple sort is one of the variants of Sorting by Selection,see Knuth's TAOCP Vol 3 Ch. 5.2.3 pp. 138 ff.Hugo Pfoertner=== === Subject: : Re: Sorting One of my quests is to find a source, not code, that explains the difference between a ripple sort and a bubble sort. I have Sedgewick's book, no mention of ripple sort in the index. No mention of ripple sort in the index of the Knuth edition I have. I have other, more interesting quests, but wish to resolve this one first.Google is your friend.Dave SeamanJudge Yohn's mistakes revealed in Mumia Abu-Jamal ruling.=== === Subject: : Re: SortingNot on this subject.Searching for Ripple sort resulted in lots of hits, many of which werereferences to stuff not in English.In any case, I want to find an authoritative reference in a book.All of the sorting books, such as Harold Lorin's, that I know about are outof print.http://www.standards.com/; See Howard Kaikow's web site.> One of my quests is to find a source, not code, that explains thedifference> between a ripple sort and a bubble sort.> I have Sedgewick's book, no mention of ripple sort in the index.> No mention of ripple sort in the index of the Knuth edition I have.> I have other, more interesting quests, but wish to resolve this onefirst. Google is your friend. -- Judge Yohn's mistakes revealed in Mumia Abu-Jamal ruling.