Locating and Classifying Critical Points of a 2D Function via Particle Swarm Optimization
Housam Binous
Locating and Classifying Critical Points of a 2D Function via Particle Swarm Optimization
Housam Binous
Housam Binous
This notebook implements a Particle Swarm Optimization (PSO) routine and uses it as to determine all the critical points — minima, maxima, and saddle points — of the two-variable test function f(x,y). Rather than solving ∇f = 0 symbolically (which can miss solutions or become intractable for more complex functions), the notebook redifines the problem as an optimization task: it builds a loss function equal to the squared norm of the gradient, ‖∇f‖², and searches for points where this loss is (near) zero. Since a zero gradient is exactly the condition for a critical point, any location the swarm converges to with a near-zero loss is a critical point of f(x,y).
The notebook is uses the psoOptimize function, supported by two helper routines: clipMatrix and clampMatrix that keep particle positions within bounds and velocities within a maximum speed. The swarm follows the standard PSO update rules — each particle is pulled toward its own best-known position and the swarm’s global best position, with an inertia weight that linearly decreases from 0.9 to 0.4 over the run (encouraging broad exploration early and fine-tuned convergence later) and acceleration coefficients of 1.42 for both the cognitive and social terms. The gradient and Hessian of f(x,y) are computed symbolically using D[f[x,y], {vars}] and Outer (applying D pairwise across variables), so the loss and its curvature information are always in accordance with the symbolic definition of f(x,y).
Because PSO is a stochastic algorithm that only finds one critical point per run — and a single run does not guarante to find every critical point in the domain — the notebook runs the swarm many times independently with different random seeds to increase the likelyhood of discovering all critical points. After each run converges, a classify function evaluates the Hessian’s eigenvalues at the resulting point: if all eigenvalues are positive, it’s labeled a minimum; if all negative, a maximum; and if the eigenvalues have mixed signs, it’s a saddle point. Results are printed in a table showing each run’s converged (x,y) location, the final loss value, and its classification.
Finally, the notebook produces a ContourPlot of f(x,y) over the search domain, with all the discovered critical points overlaid as colored dots — blue for minima, black for maxima, and dark green for saddle points — giving a clear visual map of where each type of critical point sits relative to the function’s peaks, valleys, and ridges.
As shown in the coutour plots, the combination of stochastic global search, automatic differentiation, and eigenvalue-based classification successfully yields all critical points of the 2D test function f(x,y).
This notebook implements a Particle Swarm Optimization (PSO) routine and uses it as to determine all the critical points — minima, maxima, and saddle points — of the two-variable test function f(x,y). Rather than solving ∇f = 0 symbolically (which can miss solutions or become intractable for more complex functions), the notebook redifines the problem as an optimization task: it builds a loss function equal to the squared norm of the gradient, ‖∇f‖², and searches for points where this loss is (near) zero. Since a zero gradient is exactly the condition for a critical point, any location the swarm converges to with a near-zero loss is a critical point of f(x,y).
The notebook is uses the psoOptimize function, supported by two helper routines: clipMatrix and clampMatrix that keep particle positions within bounds and velocities within a maximum speed. The swarm follows the standard PSO update rules — each particle is pulled toward its own best-known position and the swarm’s global best position, with an inertia weight that linearly decreases from 0.9 to 0.4 over the run (encouraging broad exploration early and fine-tuned convergence later) and acceleration coefficients of 1.42 for both the cognitive and social terms. The gradient and Hessian of f(x,y) are computed symbolically using D[f[x,y], {vars}] and Outer (applying D pairwise across variables), so the loss and its curvature information are always in accordance with the symbolic definition of f(x,y).
Because PSO is a stochastic algorithm that only finds one critical point per run — and a single run does not guarante to find every critical point in the domain — the notebook runs the swarm many times independently with different random seeds to increase the likelyhood of discovering all critical points. After each run converges, a classify function evaluates the Hessian’s eigenvalues at the resulting point: if all eigenvalues are positive, it’s labeled a minimum; if all negative, a maximum; and if the eigenvalues have mixed signs, it’s a saddle point. Results are printed in a table showing each run’s converged (x,y) location, the final loss value, and its classification.
Finally, the notebook produces a ContourPlot of f(x,y) over the search domain, with all the discovered critical points overlaid as colored dots — blue for minima, black for maxima, and dark green for saddle points — giving a clear visual map of where each type of critical point sits relative to the function’s peaks, valleys, and ridges.
As shown in the coutour plots, the combination of stochastic global search, automatic differentiation, and eigenvalue-based classification successfully yields all critical points of the 2D test function f(x,y).
In[]:=
(*ParticleSwarmOptimizer*)clipMatrix[mat_,lbv_,ubv_]:=Transpose[MapThread[Clip[#1,{#2,#3}]&,{Transpose[mat],lbv,ubv}]];clampMatrix[mat_,vmaxv_]:=Transpose[MapThread[Clip[#1,{-#2,#2}]&,{Transpose[mat],vmaxv}]];
In[]:=
psoOptimize[lossFunc_,dim_,lb_,ub_,nParticles_:30,nIter_:100,wMax_:0.9,wMin_:0.4,c1_:1.42,c2_:1.42,vClampFrac_:0.2,seed_:1]:=Module[{particles,velocities,pbest,pbestScore,gbest,gbestScore,convergence={},r1,r2,scores,gbestIdx,vmax,wIter},SeedRandom[seed];vmax=vClampFrac*(ub-lb);particles=Transpose@Table[RandomReal[{lb[[i]],ub[[i]]},nParticles],{i,dim}];velocities=ConstantArray[0.,{nParticles,dim}];pbest=particles;pbestScore=lossFunc/@particles;gbestIdx=First@Ordering[pbestScore,1];gbestScore=pbestScore[[gbestIdx]];gbest=pbest[[gbestIdx]];Do[wIter=wMax-(wMax-wMin)*(j-1)/(nIter-1);r1=RandomReal[1,{nParticles,dim}];r2=RandomReal[1,{nParticles,dim}];velocities=wIter*velocities+c1*r1*(pbest-particles)+c2*r2*(ConstantArray[gbest,nParticles]-particles);velocities=clampMatrix[velocities,vmax];particles+=velocities;particles=clipMatrix[particles,lb,ub];scores=lossFunc/@particles;pbest=MapThread[If[#3<#2,#1,#4]&,{particles,pbestScore,scores,pbest}];pbestScore=MapThread[Min,{pbestScore,scores}];gbestIdx=First@Ordering[pbestScore,1];If[pbestScore[[gbestIdx]]<gbestScore,gbestScore=pbestScore[[gbestIdx]];gbest=pbest[[gbestIdx]];];AppendTo[convergence,gbestScore],{j,nIter}];{gbest,gbestScore}];
In[]:=
(*ClassifyeachconvergedpointviatheHessianeigenvalues*)classify[{x_,y_}]:=Module[{eig},eig=Eigenvalues[hessf[x,y]];Which[Max[Abs[eig]]>1*^-4&&Min[eig]>0,"minimum",Max[Abs[eig]]>1*^-4&&Max[eig]<0,"maximum",Min[eig]<0&&Max[eig]>0,"saddle"]];
Test funtion 1
In[]:=
(*FindingallcriticalpointsviaPSObyminimizing||gradf||^2Testfunction:f(x,y)=Sin[2Pix]*Cos[3Piy]*)ClearAll[f,loss,gradf,hessf];f[x_,y_]:=Sin[2Pix]*Cos[3Piy];
In[]:=
vars={x,y};gradf[x_,y_]:=Evaluate@D[f[x,y],{vars}](*Gradient:{∂f/∂x,∂f/∂y}*)hessf[x_,y_]:=Evaluate@Outer[D[f[x,y],{#1,1},{#2,1}]&,vars,vars](*HessianusingOuter*)
In[]:=
gradf[x,y]
Out[]=
{2πCos[2πx]Cos[3πy],-3πSin[2πx]Sin[3πy]}
In[]:=
hessf[x,y]//MatrixForm
Out[]//MatrixForm=
-4 2 π | -6 2 π |
-6 2 π | -9 2 π |
In[]:=
(*ObjectivePSOwillminimize:squarednormofthegradient*)loss[{x_,y_}]:=Total[gradf[x,y]^2];
In[]:=
(*Runseveralindependentswarms(differentrandomseeds)sowe'relikelytolandonallcriticalpoints*)results=Table[psoOptimize[loss,2,{0,0},{1,1},30,100,0.9,0.4,1.42,1.42,0.2,seed],{seed,1,100}];
In[]:=
Print["Run gbest (x,y) ||grad||^2 classification"];Do[With[{pt=results[[k,1]],val=results[[k,2]]},Print[k," ",NumberForm[pt,{6,4}]," ",ScientificForm[val,3]," ",classify[pt]]],{k,Length[results]}];
Run gbest (x,y) ||grad||^2 classification
1 {0.5000,0.1667} 1.83× saddle
-16
10
2 {1.0000,0.5000} 4.51× saddle
-19
10
3 {1.0000,0.8333} 6.69× saddle
-20
10
4 {0.2500,1.0000} 8.28× minimum
-22
10
5 {0.0000,0.5000} 2.47× saddle
-20
10
6 {0.2500,0.6667} 6.51× maximum
-15
10
7 {0.7500,0.3333} 3.44× maximum
-14
10
8 {1.0000,0.5000} 1.56× saddle
-20
10
9 {0.2500,0.0000} 6.81× maximum
-19
10
10 {0.5000,0.5000} 1.47× saddle
-17
10
11 {0.5000,0.5000} 7.51× saddle
-17
10
12 {0.7500,0.3333} 6.05× maximum
-17
10
13 {0.5000,0.5000} 3.25× saddle
-16
10
14 {0.7500,0.3333} 7.06× maximum
-16
10
15 {1.0000,0.1667} 2.19× saddle
-16
10
16 {0.0000,0.5000} 3.8× saddle
-22
10
17 {0.5000,0.5000} 7.64× saddle
-15
10
18 {0.7500,0.0000} 5.26× minimum
-20
10
19 {0.5000,0.5000} 1.11× saddle
-16
10
20 {1.0000,0.8333} 7.17× saddle
-19
10
Test funtion 2