One of the typical aims of statistical shape analysis, in addition to deriving an estimate of mean shape, is to get an estimate of shape variability. This aim is achived through employing the principal component analysis. Because the principal component analysis is limited to data on Euclidean space, this method cannot be applied for the shape data which are inherently non-Euclidean data. In this situation, the principal geodesic analysis or its linear approximation can be used as a generalization of the principal component analysis in non-Euclidean space. Because the main root of this method is the gradient descent algorithm, revealing some of its main defects, a new algorithm is proposed in this paper which leads to a robust estimate of mean shape and also preserves the geometrical structure of shape. Then, providing some theoretical aspects of principal geodesic analysis, its application is evaluated in a simulation study and in real data.