This paper presents a simple numerical procedure based upon the projected gradient descent (PGD) and finite element method (FEM) for the shape optimization of laterally restrained columns to attain the maximum elastic buckling load under the specified volumetric constraint. The analysis of the buckling load is achieved via the formulation based on Euler-Bernoulli beam theory, the discretization by the standard finite element technique, and the determination of the least eigenvalue and the corresponding eigenvector via the power method with Rayleigh quotient. In the optimization, the profile of the cross-sectional area of the column is represented by piecewise polynomial interpolation functions. The gradient information and the projection operator required in PGD iterations are obtained explicitly in a closed form. A selected set of results is reported to demonstrate not only the good convergence behavior and accuracy of numerical solutions, but also the capability of the proposed technique to attain the optimal shape of columns for various scenarios.