Hessian Matrix and Its Role in Convex Optimization

mathematical statistics
optimization
R
The Hessian matrix is the multivariable generalization of the second derivative. It captures curvature information in all directions simultaneously.
Author

Yupng L.

Published

December 4, 2025

The Hessian Matrix and Its Role in Convex Optimization

Table of Contents

  1. Introduction: From Calculus to Optimization
  2. Definition of the Hessian Matrix
  3. Examples of Hessian Matrices
  4. Taylor Expansion and the Hessian
  5. The Hessian and Convexity
  6. Optimization: Finding Minima
  7. The Hessian in Optimization Algorithms
  8. Applications in Statistics and Machine Learning
  9. Computational Considerations
  10. Summary

1. Introduction: From Calculus to Optimization

In single-variable calculus, the second derivative \(f''(x)\) tells us about the curvature of a function:

  • \(f''(x) > 0\): function curves upward (convex locally)
  • \(f''(x) < 0\): function curves downward (concave locally)
  • \(f''(x) = 0\): possible inflection point

The Hessian matrix is the multivariable generalization of the second derivative. It captures curvature information in all directions simultaneously.

Why Does This Matter?

Single Variable Multivariable
\(f''(x) > 0\) at critical point → local minimum \(\mathbf{H}_f(\mathbf{x}) \succ 0\) at critical point → local minimum
\(f''(x) > 0\) everywhere → convex function \(\mathbf{H}_f(\mathbf{x}) \succeq 0\) everywhere → convex function
\(f''(x) < 0\) at critical point → local maximum \(\mathbf{H}_f(\mathbf{x}) \prec 0\) at critical point → local maximum

2. Definition of the Hessian Matrix

Definition (Hessian Matrix)

Let \(f: \mathbb{R}^n \to \mathbb{R}\) be a twice continuously differentiable function. The Hessian matrix of \(f\) at point \(\mathbf{x} \in \mathbb{R}^n\) is the \(n \times n\) matrix of second-order partial derivatives:

\[ \mathbf{H}_f(\mathbf{x}) = \nabla^2 f(\mathbf{x}) = \begin{pmatrix} \dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_1 \partial x_n} \\[12pt] \dfrac{\partial^2 f}{\partial x_2 \partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \cdots & \dfrac{\partial^2 f}{\partial x_2 \partial x_n} \\[12pt] \vdots & \vdots & \ddots & \vdots \\[12pt] \dfrac{\partial^2 f}{\partial x_n \partial x_1} & \dfrac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_n^2} \end{pmatrix} \]

In index notation:

\[ [\mathbf{H}_f(\mathbf{x})]_{ij} = \frac{\partial^2 f}{\partial x_i \partial x_j} \bigg|_{\mathbf{x}} \]

Theorem (Symmetry of the Hessian)

If \(f\) has continuous second partial derivatives (i.e., \(f \in C^2\)), then by Schwarz’s theorem (also called Clairaut’s theorem):

\[ \frac{\partial^2 f}{\partial x_i \partial x_j} = \frac{\partial^2 f}{\partial x_j \partial x_i} \]

Therefore, \(\mathbf{H}_f(\mathbf{x})\) is a symmetric matrix.

Definition (Gradient Vector)

For completeness, the gradient of \(f\) is the vector of first partial derivatives:

\[ \nabla f(\mathbf{x}) = \begin{pmatrix} \dfrac{\partial f}{\partial x_1} \\[10pt] \dfrac{\partial f}{\partial x_2} \\[10pt] \vdots \\[10pt] \dfrac{\partial f}{\partial x_n} \end{pmatrix} \]

The Hessian can be viewed as the Jacobian of the gradient:

\[ \mathbf{H}_f(\mathbf{x}) = \mathbf{J}(\nabla f)(\mathbf{x}) \]

3. Examples of Hessian Matrices

Example 1: Quadratic Function

Consider the quadratic function:

\[ f(\mathbf{x}) = \frac{1}{2}\mathbf{x}^\top \mathbf{A} \mathbf{x} - \mathbf{b}^\top \mathbf{x} + c \]

where \(\mathbf{A}\) is an \(n \times n\) symmetric matrix, \(\mathbf{b} \in \mathbb{R}^n\), and \(c \in \mathbb{R}\).

Gradient:

\[ \nabla f(\mathbf{x}) = \mathbf{A}\mathbf{x} - \mathbf{b} \]

Hessian:

\[ \mathbf{H}_f(\mathbf{x}) = \mathbf{A} \]

Note: The Hessian is constant (independent of \(\mathbf{x}\)) for quadratic functions. This is why quadratic functions are the “simplest” case for optimization analysis.

Example 2: Two-Variable Function

Let \(f(x_1, x_2) = x_1^3 + 2x_1^2 x_2 - x_2^2 + 3x_1\).

First derivatives (gradient):

\[ \frac{\partial f}{\partial x_1} = 3x_1^2 + 4x_1 x_2 + 3 \]

\[ \frac{\partial f}{\partial x_2} = 2x_1^2 - 2x_2 \]

Second derivatives:

\[ \frac{\partial^2 f}{\partial x_1^2} = 6x_1 + 4x_2 \]

\[ \frac{\partial^2 f}{\partial x_1 \partial x_2} = 4x_1 \]

\[ \frac{\partial^2 f}{\partial x_2^2} = -2 \]

Hessian matrix:

\[ \mathbf{H}_f(x_1, x_2) = \begin{pmatrix} 6x_1 + 4x_2 & 4x_1 \\ 4x_1 & -2 \end{pmatrix} \]

Note: This Hessian depends on \(\mathbf{x}\), so the curvature changes at different points.

Example 3: Log-Sum-Exp Function

The log-sum-exp function is fundamental in machine learning (softmax, logistic regression):

\[ f(\mathbf{x}) = \log\left(\sum_{i=1}^n e^{x_i}\right) \]

Let \(z_i = e^{x_i}\) and \(Z = \sum_i z_i\). Then:

Gradient:

\[ \frac{\partial f}{\partial x_i} = \frac{z_i}{Z} = p_i \quad \text{(softmax probabilities)} \]

Hessian entries:

\[ \frac{\partial^2 f}{\partial x_i \partial x_j} = \begin{cases} p_i(1 - p_i) & \text{if } i = j \\[8pt] -p_i p_j & \text{if } i \neq j \end{cases} \]

In matrix form, letting \(\mathbf{p} = (p_1, \ldots, p_n)^\top\):

\[ \mathbf{H}_f(\mathbf{x}) = \text{diag}(\mathbf{p}) - \mathbf{p}\mathbf{p}^\top \]

This Hessian is positive semi-definite, proving that log-sum-exp is convex!

Example 4: Least Squares Loss

For linear regression with loss function:

\[ f(\boldsymbol{\beta}) = \frac{1}{2}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 = \frac{1}{2}(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \]

Gradient:

\[ \nabla f(\boldsymbol{\beta}) = -\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) = \mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} - \mathbf{X}^\top\mathbf{y} \]

Hessian:

\[ \mathbf{H}_f(\boldsymbol{\beta}) = \mathbf{X}^\top\mathbf{X} \]

Since \(\mathbf{X}^\top\mathbf{X}\) is always PSD (and PD when \(\mathbf{X}\) has full column rank), the least squares loss is convex!

4. Taylor Expansion and the Hessian

The fundamental connection between the Hessian and optimization comes from the Taylor expansion.

Theorem (Second-Order Taylor Expansion)

For \(f \in C^2(\mathbb{R}^n)\), around point \(\mathbf{x}\):

\[ \boxed{f(\mathbf{x} + \mathbf{d}) = f(\mathbf{x}) + \nabla f(\mathbf{x})^\top \mathbf{d} + \frac{1}{2}\mathbf{d}^\top \mathbf{H}_f(\mathbf{x}) \mathbf{d} + O(\|\mathbf{d}\|^3)} \]

where \(\mathbf{d} \in \mathbb{R}^n\) is a displacement vector.

Interpretation of Each Term

Term Name Meaning
\(f(\mathbf{x})\) Function value Starting point value
\(\nabla f(\mathbf{x})^\top \mathbf{d}\) Linear term First-order (linear) approximation
\(\frac{1}{2}\mathbf{d}^\top \mathbf{H}_f(\mathbf{x}) \mathbf{d}\) Quadratic term Second-order curvature correction
\(O(\|\mathbf{d}\|^3)\) Remainder Higher-order terms (negligible for small \(\mathbf{d}\))

Directional Curvature

The term \(\mathbf{d}^\top \mathbf{H}_f(\mathbf{x}) \mathbf{d}\) captures how the function curves in direction \(\mathbf{d}\).

Definition (Directional Second Derivative):

For a unit vector \(\mathbf{u}\) (\(\|\mathbf{u}\| = 1\)), the second directional derivative of \(f\) at \(\mathbf{x}\) in direction \(\mathbf{u}\) is:

\[ \frac{d^2}{dt^2} f(\mathbf{x} + t\mathbf{u}) \bigg|_{t=0} = \mathbf{u}^\top \mathbf{H}_f(\mathbf{x}) \mathbf{u} \]

This is the curvature of \(f\) along direction \(\mathbf{u}\).

Connection to Eigenvalues

If \(\mathbf{H}_f(\mathbf{x})\) has eigenvalues \(\lambda_1 \leq \lambda_2 \leq \cdots \leq \lambda_n\):

  • Maximum curvature = \(\lambda_n\) (along eigenvector \(\mathbf{v}_n\))
  • Minimum curvature = \(\lambda_1\) (along eigenvector \(\mathbf{v}_1\))

For any unit vector \(\mathbf{u}\):

\[ \lambda_1 \leq \mathbf{u}^\top \mathbf{H}_f(\mathbf{x}) \mathbf{u} \leq \lambda_n \]

5. The Hessian and Convexity

5.1 Convex Functions: Definition

Definition (Convex Function):

A function \(f: \mathbb{R}^n \to \mathbb{R}\) is convex if for all \(\mathbf{x}, \mathbf{y} \in \mathbb{R}^n\) and \(\lambda \in [0, 1]\):

\[ f(\lambda \mathbf{x} + (1-\lambda)\mathbf{y}) \leq \lambda f(\mathbf{x}) + (1-\lambda) f(\mathbf{y}) \]

Geometric interpretation: The line segment connecting any two points on the graph lies above (or on) the graph.

Definition (Strictly Convex Function):

\(f\) is strictly convex if the inequality is strict for \(\lambda \in (0, 1)\) and \(\mathbf{x} \neq \mathbf{y}\).

5.2 The Fundamental Theorem: Hessian Characterization of Convexity

Theorem (Hessian Test for Convexity):

Let \(f: \mathbb{R}^n \to \mathbb{R}\) be twice continuously differentiable on a convex set \(\mathcal{C}\). Then:

(a) \(f\) is convex on \(\mathcal{C}\) if and only if:

\[ \boxed{\mathbf{H}_f(\mathbf{x}) \succeq \mathbf{0} \quad \forall \mathbf{x} \in \mathcal{C}} \]

(The Hessian is positive semi-definite everywhere)

(b) \(f\) is strictly convex on \(\mathcal{C}\) if:

\[ \boxed{\mathbf{H}_f(\mathbf{x}) \succ \mathbf{0} \quad \forall \mathbf{x} \in \mathcal{C}} \]

(The Hessian is positive definite everywhere)

(c) \(f\) is concave on \(\mathcal{C}\) if and only if \(\mathbf{H}_f(\mathbf{x}) \preceq \mathbf{0}\) (negative semi-definite)

5.3 Proof Sketch

(\(\Rightarrow\)) If \(f\) is convex, then \(\mathbf{H}_f(\mathbf{x}) \succeq 0\):

For any \(\mathbf{x}\) and direction \(\mathbf{d}\), define \(g(t) = f(\mathbf{x} + t\mathbf{d})\).

Since \(f\) is convex, \(g\) is a convex function of one variable. Thus:

\[ g''(0) = \mathbf{d}^\top \mathbf{H}_f(\mathbf{x}) \mathbf{d} \geq 0 \]

Since \(\mathbf{d}\) was arbitrary, \(\mathbf{H}_f(\mathbf{x}) \succeq \mathbf{0}\).

(\(\Leftarrow\)) If \(\mathbf{H}_f(\mathbf{x}) \succeq 0\) everywhere, then \(f\) is convex:

By Taylor’s theorem with integral remainder:

\[ f(\mathbf{y}) = f(\mathbf{x}) + \nabla f(\mathbf{x})^\top (\mathbf{y} - \mathbf{x}) + \int_0^1 (1-t)(\mathbf{y}-\mathbf{x})^\top \mathbf{H}_f(\mathbf{x} + t(\mathbf{y}-\mathbf{x}))(\mathbf{y}-\mathbf{x})\, dt \]

Since the Hessian is PSD, the integral is non-negative:

\[ f(\mathbf{y}) \geq f(\mathbf{x}) + \nabla f(\mathbf{x})^\top (\mathbf{y} - \mathbf{x}) \]

This is the first-order condition for convexity, which implies convexity.

5.4 Summary Table: Hessian Properties and Function Behavior

Hessian Property Notation Function Property Geometric Shape
Positive Definite (PD) \(\mathbf{H} \succ 0\) Strictly convex Bowl-shaped, curves up in ALL directions
Positive Semi-Definite (PSD) \(\mathbf{H} \succeq 0\) Convex Bowl or flat, never curves down
Negative Definite (ND) \(\mathbf{H} \prec 0\) Strictly concave Dome-shaped, curves down in ALL directions
Negative Semi-Definite (NSD) \(\mathbf{H} \preceq 0\) Concave Dome or flat, never curves up
Indefinite Mixed eigenvalues Neither convex nor concave Saddle point shape

6. Optimization: Finding Minima

6.1 First-Order Necessary Condition

Theorem: If \(\mathbf{x}^*\) is a local minimum of \(f\) and \(f\) is differentiable at \(\mathbf{x}^*\), then:

\[ \nabla f(\mathbf{x}^*) = \mathbf{0} \]

Points satisfying this are called stationary points or critical points.

6.2 Second-Order Conditions

Theorem (Second-Order Necessary Condition):

If \(\mathbf{x}^*\) is a local minimum and \(f \in C^2\), then:

\[ \nabla f(\mathbf{x}^*) = \mathbf{0} \quad \text{and} \quad \mathbf{H}_f(\mathbf{x}^*) \succeq \mathbf{0} \]

Theorem (Second-Order Sufficient Condition):

If at point \(\mathbf{x}^*\):

\[ \nabla f(\mathbf{x}^*) = \mathbf{0} \quad \text{and} \quad \mathbf{H}_f(\mathbf{x}^*) \succ \mathbf{0} \]

then \(\mathbf{x}^*\) is a strict local minimum.

6.3 Classification of Critical Points

At a critical point \(\mathbf{x}^*\) where \(\nabla f(\mathbf{x}^*) = \mathbf{0}\):

Hessian at \(\mathbf{x}^*\) Eigenvalues Classification
\(\mathbf{H} \succ 0\) (PD) All \(\lambda_i > 0\) Strict local minimum
\(\mathbf{H} \prec 0\) (ND) All \(\lambda_i < 0\) Strict local maximum
Indefinite Some \(\lambda_i > 0\), some \(\lambda_i < 0\) Saddle point
\(\mathbf{H} \succeq 0\) or \(\mathbf{H} \preceq 0\) Some \(\lambda_i = 0\) Inconclusive (need higher-order info)

6.4 Global vs Local Minima

For general functions: A local minimum may not be global.

For convex functions: This is where convexity shines!

Theorem (Global Optimality for Convex Functions):

If \(f\) is convex and \(\mathbf{x}^*\) is a local minimum, then \(\mathbf{x}^*\) is a global minimum.

Moreover, if \(f\) is strictly convex, the global minimum (if it exists) is unique.

This is the key reason why convexity matters in optimization!

7. The Hessian in Optimization Algorithms

7.1 Newton’s Method

Newton’s method uses the Hessian to find the minimum by iteratively solving:

\[ \mathbf{x}_{k+1} = \mathbf{x}_k - [\mathbf{H}_f(\mathbf{x}_k)]^{-1} \nabla f(\mathbf{x}_k) \]

Derivation: At each step, minimize the second-order Taylor approximation:

\[ f(\mathbf{x}_k + \mathbf{d}) \approx f(\mathbf{x}_k) + \nabla f(\mathbf{x}_k)^\top \mathbf{d} + \frac{1}{2}\mathbf{d}^\top \mathbf{H}_f(\mathbf{x}_k) \mathbf{d} \]

Taking the derivative with respect to \(\mathbf{d}\) and setting to zero:

\[ \nabla f(\mathbf{x}_k) + \mathbf{H}_f(\mathbf{x}_k) \mathbf{d} = \mathbf{0} \]

\[ \mathbf{d}^* = -[\mathbf{H}_f(\mathbf{x}_k)]^{-1} \nabla f(\mathbf{x}_k) \]

Properties of Newton’s Method:

Aspect Description
Convergence Quadratic (very fast) near minimum
Requirement Hessian must be PD for descent direction
Cost \(O(n^3)\) per iteration (matrix inversion)
Storage \(O(n^2)\) for Hessian

7.2 Comparison of Optimization Methods

Method Uses Gradient? Uses Hessian? Convergence Rate Cost per Iteration
Gradient Descent Yes No Linear \(O(n)\)
Newton’s Method Yes Yes (exact) Quadratic \(O(n^3)\)
Quasi-Newton (BFGS) Yes Yes (approximate) Superlinear \(O(n^2)\)
Conjugate Gradient Yes No (uses Hessian-vector products) Linear to Superlinear \(O(n)\)

7.3 Quasi-Newton Methods

Since computing and inverting the Hessian is expensive, quasi-Newton methods build an approximation \(\mathbf{B}_k \approx \mathbf{H}_f(\mathbf{x}_k)\) iteratively.

BFGS Update Formula:

\[ \mathbf{B}_{k+1} = \mathbf{B}_k - \frac{\mathbf{B}_k \mathbf{s}_k \mathbf{s}_k^\top \mathbf{B}_k}{\mathbf{s}_k^\top \mathbf{B}_k \mathbf{s}_k} + \frac{\mathbf{y}_k \mathbf{y}_k^\top}{\mathbf{y}_k^\top \mathbf{s}_k} \]

where: - \(\mathbf{s}_k = \mathbf{x}_{k+1} - \mathbf{x}_k\) - \(\mathbf{y}_k = \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_k)\)

Key property: If \(f\) is convex (so the true Hessian is PSD), and certain conditions hold, the BFGS approximation remains PD, ensuring descent directions.

7.4 When the Hessian is Not PD

If \(\mathbf{H}_f(\mathbf{x}_k)\) is not positive definite (e.g., near a saddle point), Newton’s direction may not be a descent direction.

Solutions:

  1. Modified Newton: Replace \(\mathbf{H}\) with \(\mathbf{H} + \lambda \mathbf{I}\) where \(\lambda > 0\) makes it PD
  2. Trust Region Methods: Restrict step size to region where quadratic model is valid
  3. Line Search: Ensure sufficient decrease along the Newton direction

8. Applications in Statistics and Machine Learning

8.1 Maximum Likelihood Estimation (MLE)

For log-likelihood \(\ell(\boldsymbol{\theta})\), we want to find:

\[ \hat{\boldsymbol{\theta}}_{\text{MLE}} = \arg\max_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}) \]

Newton-Raphson for MLE:

\[ \boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k - [\mathbf{H}_\ell(\boldsymbol{\theta}_k)]^{-1} \nabla \ell(\boldsymbol{\theta}_k) \]

Fisher Scoring: Replace observed Hessian with expected Hessian (Fisher Information):

\[ \mathbf{I}(\boldsymbol{\theta}) = -\mathbb{E}[\mathbf{H}_\ell(\boldsymbol{\theta})] = \mathbb{E}\left[\nabla \ell(\boldsymbol{\theta}) \nabla \ell(\boldsymbol{\theta})^\top\right] \]

8.2 Logistic Regression

For logistic regression with log-likelihood:

\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ y_i \mathbf{x}_i^\top \boldsymbol{\beta} - \log(1 + e^{\mathbf{x}_i^\top \boldsymbol{\beta}}) \right] \]

Gradient:

\[ \nabla \ell(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i - p_i) \mathbf{x}_i = \mathbf{X}^\top (\mathbf{y} - \mathbf{p}) \]

where \(p_i = \frac{1}{1 + e^{-\mathbf{x}_i^\top \boldsymbol{\beta}}}\)

Hessian:

\[ \mathbf{H}_\ell(\boldsymbol{\beta}) = -\sum_{i=1}^n p_i(1-p_i) \mathbf{x}_i \mathbf{x}_i^\top = -\mathbf{X}^\top \mathbf{W} \mathbf{X} \]

where \(\mathbf{W} = \text{diag}(p_1(1-p_1), \ldots, p_n(1-p_n))\)

Key insight: Since \(0 < p_i < 1\), we have \(p_i(1-p_i) > 0\), so \(\mathbf{W} \succ 0\). Thus \(\mathbf{H}_\ell \prec 0\), meaning the log-likelihood is concave and has a unique maximum!

8.3 Neural Networks and Non-Convexity

Neural network loss functions are typically non-convex:

  • The Hessian can be indefinite (saddle points are common)
  • Multiple local minima exist
  • Gradient descent can get stuck

Role of the Hessian in deep learning:

  1. Understanding loss landscape: Eigenvalues of Hessian reveal curvature structure
  2. Second-order methods: Natural gradient, K-FAC approximate Hessian efficiently
  3. Saddle point detection: Negative eigenvalues indicate escape directions

8.4 Standard Errors and Confidence Intervals

In MLE, the asymptotic covariance of \(\hat{\boldsymbol{\theta}}\) is:

\[ \text{Cov}(\hat{\boldsymbol{\theta}}) \approx [-\mathbf{H}_\ell(\hat{\boldsymbol{\theta}})]^{-1} = \mathbf{I}(\hat{\boldsymbol{\theta}})^{-1} \]

This requires \(-\mathbf{H}_\ell\) to be positive definite (i.e., \(\ell\) is concave at the MLE).

Standard errors: \(\text{SE}(\hat{\theta}_j) = \sqrt{[\mathbf{I}(\hat{\boldsymbol{\theta}})^{-1}]_{jj}}\)

9. Computational Considerations

9.1 Checking Positive Definiteness

Methods to check if \(\mathbf{H}\) is PD:

Method Complexity Notes
Compute all eigenvalues \(O(n^3)\) Check if all \(\lambda_i > 0\)
Cholesky decomposition \(O(n^3/3)\) Succeeds iff PD
Leading principal minors \(O(n^4)\) Check if all positive
Sylvester’s criterion \(O(n^4)\) Compute \(n\) determinants

In practice: Try Cholesky decomposition. If it fails, the matrix is not PD.

9.2 Numerical Issues

Problem: Computed Hessian may not be exactly symmetric or PD due to: - Floating-point errors - Finite difference approximations - Numerical instability

Solutions:

  1. Symmetrize: \(\mathbf{H} \leftarrow \frac{1}{2}(\mathbf{H} + \mathbf{H}^\top)\)

  2. Regularize: \(\mathbf{H} \leftarrow \mathbf{H} + \epsilon \mathbf{I}\) for small \(\epsilon > 0\)

  3. Project to PSD cone: Use eigendecomposition, clip negative eigenvalues to zero (or small positive value)

9.3 Hessian-Free Methods

For large-scale problems, storing the \(n \times n\) H

Reference

Content generated by Claude (Claude Opus 4.5), an AI assistant developed by Anthropic. Conversation date: [Your Date]. Retrieved from: https://claude.ai/

Appendix - R Code for Illustration

# =============================================================================
# HESSIAN MATRIX AND CONVEX OPTIMIZATION IN R
# =============================================================================

library(ggplot2)
library(numDeriv)  # For numerical derivatives

# -----------------------------------------------------------------------------
# 1. COMPUTING THE HESSIAN
# -----------------------------------------------------------------------------

# Example function: f(x,y) = x^2 + 2y^2 + xy
f <- function(x) x[1]^2 + 2*x[2]^2 + x[1]*x[2]

# Analytical gradient
grad_f <- function(x) c(2*x[1] + x[2], 4*x[2] + x[1])

# Analytical Hessian (constant for quadratic)
hess_f <- function(x) matrix(c(2, 1, 1, 4), nrow = 2)

# Numerical Hessian (for verification)
x0 <- c(1, 1)
H_numerical <- hessian(f, x0)
H_analytical <- hess_f(x0)

cat("Analytical Hessian:\n"); print(H_analytical)
Analytical Hessian:
     [,1] [,2]
[1,]    2    1
[2,]    1    4
cat("\nNumerical Hessian:\n"); print(round(H_numerical, 6))

Numerical Hessian:
     [,1] [,2]
[1,]    2    1
[2,]    1    4
# -----------------------------------------------------------------------------
# 2. CHECK POSITIVE DEFINITENESS
# -----------------------------------------------------------------------------

check_definiteness <- function(H) {
  eig <- eigen(H, symmetric = TRUE)$values
  if (all(eig > 1e-10)) return("Positive Definite (Convex)")
  if (all(eig >= -1e-10)) return("Positive Semi-Definite")
  if (all(eig < -1e-10)) return("Negative Definite (Concave)")
  return("Indefinite (Saddle Point)")
}

cat("\nDefiniteness:", check_definiteness(H_analytical))

Definiteness: Positive Definite (Convex)
cat("\nEigenvalues:", eigen(H_analytical)$values, "\n")

Eigenvalues: 4.414214 1.585786 
# -----------------------------------------------------------------------------
# 3. VISUALIZE CONVEX FUNCTION
# -----------------------------------------------------------------------------

# Create contour plot
grid <- expand.grid(x = seq(-3, 3, 0.1), y = seq(-3, 3, 0.1))
grid$z <- apply(grid, 1, f)

ggplot(grid, aes(x, y, z = z)) +
  geom_contour_filled(bins = 20) +
  geom_point(aes(x = 0, y = 0), color = "red", size = 3) +
  annotate("text", x = 0.3, y = 0.3, label = "Minimum", color = "red") +
  labs(title = "Convex Function (PD Hessian)", 
       subtitle = "Unique global minimum at origin") +
  theme_minimal() + coord_fixed()

# -----------------------------------------------------------------------------
# 4. NEWTON'S METHOD
# -----------------------------------------------------------------------------

newton_method <- function(f, grad, hess, x0, tol = 1e-8, max_iter = 100) {
  x <- x0
  history <- list(x)
  
  for (i in 1:max_iter) {
    g <- grad(x)
    H <- hess(x)
    
    # Newton step: x_new = x - H^{-1} * g
    d <- solve(H, -g)
    x <- x + d
    history[[i + 1]] <- x
    
    if (sqrt(sum(g^2)) < tol) break
  }
  
  list(minimum = x, iterations = i, path = do.call(rbind, history))
}

result <- newton_method(f, grad_f, hess_f, x0 = c(2, 2))
cat("\nNewton's Method Result:\n")

Newton's Method Result:
cat("Minimum at:", result$minimum, "\n")
Minimum at: 0 0 
cat("Iterations:", result$iterations, "\n")
Iterations: 2 
# -----------------------------------------------------------------------------
# 5. COMPARE: CONVEX VS SADDLE POINT
# -----------------------------------------------------------------------------

# Saddle function: f(x,y) = x^2 - y^2
f_saddle <- function(x) x[1]^2 - x[2]^2
H_saddle <- matrix(c(2, 0, 0, -2), nrow = 2)

cat("\nSaddle Point Hessian:\n"); print(H_saddle)

Saddle Point Hessian:
     [,1] [,2]
[1,]    2    0
[2,]    0   -2
cat("Eigenvalues:", eigen(H_saddle)$values, "\n")
Eigenvalues: 2 -2 
cat("Definiteness:", check_definiteness(H_saddle), "\n")
Definiteness: Indefinite (Saddle Point) 
# Visualize saddle
grid$z_saddle <- apply(grid[, 1:2], 1, f_saddle)
ggplot(grid, aes(x, y, z = z_saddle)) +
  geom_contour_filled(bins = 20) +
  geom_point(aes(x = 0, y = 0), color = "red", size = 3) +
  labs(title = "Saddle Point (Indefinite Hessian)",
       subtitle = "Critical point is NOT a minimum") +
  theme_minimal() + coord_fixed()

# -----------------------------------------------------------------------------
# 6. LOGISTIC REGRESSION EXAMPLE
# -----------------------------------------------------------------------------

# Simulate data
set.seed(42)
n <- 200
X <- cbind(1, rnorm(n), rnorm(n))
beta_true <- c(-0.5, 1.5, -1)
p_true <- 1 / (1 + exp(-X %*% beta_true))
y <- rbinom(n, 1, p_true)

# Negative log-likelihood (to minimize)
nll <- function(beta) {
  eta <- X %*% beta
  -sum(y * eta - log(1 + exp(eta)))
}

# Gradient
grad_nll <- function(beta) {
  p <- 1 / (1 + exp(-X %*% beta))
  -t(X) %*% (y - p)
}

# Hessian
hess_nll <- function(beta) {
  p <- 1 / (1 + exp(-X %*% beta))
  W <- diag(as.vector(p * (1 - p)))
  t(X) %*% W %*% X  # Always PSD!
}

# Fit using Newton's method
beta_hat <- newton_method(nll, grad_nll, hess_nll, x0 = c(0, 0, 0))$minimum
Warning in (function (..., deparse.level = 1) : number of columns of result is
not a multiple of vector length (arg 1)
cat("\nLogistic Regression Coefficients:\n")

Logistic Regression Coefficients:
cat("True:", beta_true, "\n")
True: -0.5 1.5 -1 
cat("Estimated:", round(beta_hat, 3), "\n")
Estimated: -0.449 1.636 -1.771 
# Verify Hessian is PD
cat("\nHessian at MLE:", check_definiteness(hess_nll(beta_hat)), "\n")

Hessian at MLE: Positive Definite (Convex) 
# Standard errors from inverse Hessian
se <- sqrt(diag(solve(hess_nll(beta_hat))))
cat("Standard Errors:", round(se, 3), "\n")
Standard Errors: 0.196 0.282 0.303