{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "e466cb3d-c036-45f7-af9f-b18b9aa22e8a",
      "metadata": {},
      "source": [
        "# Optimizers and Gradients\n",
        "\n",
        "Many quantum algorithms require the optimization of quantum circuit parameters with respect to an expectation value. CUDA-Q provides a comprehensive suite of optimization tools for hybrid quantum-classical algorithms like VQE (Variational Quantum Eigensolver).\n",
        "\n",
        "This notebook will demonstrate:\n",
        "\n",
        "1. **Built-in CUDA-Q Optimizers**: Adam, SGD, SPSA, COBYLA, NelderMead, LBFGS, and GradientDescent\n",
        "2. **Optimizer Parameters**: Detailed configuration options with defaults and tuning guidance\n",
        "3. **Gradient Strategies**: CentralDifference, ForwardDifference, and ParameterShift\n",
        "4. **Third-Party Optimizers**: Integration with SciPy\n",
        "5. **Parallel Parameter Shift**: Multi-GPU gradient computation\n",
        "\n",
        "## CUDA-Q Optimizer Overview\n",
        "\n",
        "CUDA-Q includes the following optimizers:\n",
        "\n",
        "### Gradient-Free Optimizers (no gradients required):\n",
        "- **COBYLA**: Constrained Optimization BY Linear Approximations\n",
        "- **NelderMead**: Simplex-based derivative-free optimizer\n",
        "- **SPSA**: Simultaneous Perturbation Stochastic Approximation (excellent for noisy functions)\n",
        "\n",
        "### Gradient-Based Optimizers (require gradients):\n",
        "- **Adam**: Adaptive Moment Estimation with momentum (recommended for most cases)\n",
        "- **SGD**: Stochastic Gradient Descent\n",
        "- **LBFGS**: Limited-memory BFGS quasi-Newton method\n",
        "- **GradientDescent**: Basic gradient descent\n",
        "\n",
        "First, let's set up the kernel and Hamiltonian that we'll use throughout the examples."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "a77a963a-6f5c-4751-93c8-b3ccbb5921f5",
      "metadata": {},
      "outputs": [],
      "source": [
        "import cudaq\n",
        "from cudaq import spin\n",
        "import numpy as np\n",
        "\n",
        "hamiltonian = 5.907 - 2.1433 * spin.x(0) * spin.x(1) - 2.1433 * spin.y(\n",
        "    0) * spin.y(1) + .21829 * spin.z(0) - 6.125 * spin.z(1)\n",
        "\n",
        "@cudaq.kernel\n",
        "def kernel(angles: list[float]):\n",
        "    qubits = cudaq.qvector(2)\n",
        "    x(qubits[0])\n",
        "    ry(angles[0], qubits[1])\n",
        "    x.ctrl(qubits[1], qubits[0])  \n",
        "\n",
        "initial_params = np.random.normal(0, np.pi, 1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6c3de68e-754a-4412-8063-78ee9646b4b0",
      "metadata": {},
      "source": [
        "## 1. Built-in CUDA-Q Optimizers and Gradients\n",
        "\n",
        "CUDA-Q provides several optimizers with configurable parameters. Let's explore the most commonly used optimizers: **Adam**, **SGD**, and **SPSA**."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "eea5d9f2-ade2-4411-88a7-b8c5505c29c3",
      "metadata": {},
      "source": [
        "### 1.1 Adam Optimizer with Parameter Configuration\n",
        "\n",
        "**Adam (Adaptive Moment Estimation)** combines momentum and adaptive learning rates for efficient optimization. It's particularly effective for problems with noisy gradients.\n",
        "\n",
        "**Configurable Parameters:**\n",
        "- `step_size` (default: 0.01): Learning rate for parameter updates\n",
        "- `beta1` (default: 0.9): Exponential decay rate for first moment (momentum)\n",
        "- `beta2` (default: 0.999): Exponential decay rate for second moment (adaptive learning)\n",
        "- `epsilon` (default: 1e-8): Small constant for numerical stability\n",
        "- `batch_size` (default: 1): Number of samples per batch\n",
        "- `f_tol` (default: 1e-4): Convergence tolerance\n",
        "- `max_iterations`: Maximum number of iterations\n",
        "- `initial_parameters`: Starting parameter values\n",
        "\n",
        "The optimizer and gradient are specified below. An objective function is defined which uses a lambda expression to evaluate the cost (a CUDA-Q `observe` expectation value). The gradient is calculated using the `compute` method."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "87387c0a-86cf-4dfb-9130-b1c3054e43d5",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Configure Adam optimizer with custom parameters\n",
        "optimizer = cudaq.optimizers.Adam()\n",
        "optimizer.step_size = 0.1                      # Learning rate\n",
        "optimizer.beta1 = 0.9                          # First moment decay\n",
        "optimizer.beta2 = 0.999                        # Second moment decay\n",
        "optimizer.epsilon = 1e-8                       # Numerical stability\n",
        "optimizer.max_iterations = 100                 # Maximum iterations\n",
        "optimizer.initial_parameters = initial_params  # Set initial parameters\n",
        "\n",
        "# Use CentralDifference gradient strategy\n",
        "gradient = cudaq.gradients.CentralDifference()\n",
        "\n",
        "def objective_function(parameter_vector: list[float],\n",
        "                       hamiltonian=hamiltonian,\n",
        "                       gradient_strategy=gradient,\n",
        "                       kernel=kernel) -> tuple[float, list[float]]:\n",
        "    \"\"\"\n",
        "    Objective function for gradient-based optimizers.\n",
        "    Returns: (cost, gradient_vector)\n",
        "    \"\"\"\n",
        "    get_result = lambda parameter_vector: cudaq.observe(kernel, hamiltonian, parameter_vector).expectation()\n",
        "\n",
        "    cost = get_result(parameter_vector)\n",
        "    \n",
        "    gradient_vector = gradient_strategy.compute(parameter_vector, get_result, cost)\n",
        "\n",
        "    return cost, gradient_vector"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b98f7129-d360-4e9c-9ea5-d4173ecb5ef1",
      "metadata": {},
      "source": [
        "Now run the optimizer to find the optimal energy and parameters. Adam will use adaptive learning rates for each parameter."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "713d7619-7b62-42b7-bb46-601ed9883e6a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "=== Adam Optimizer Results ===\n",
            "Minimized <H> = -1.744713\n",
            "Optimal parameters: [-5.721116]\n"
          ]
        }
      ],
      "source": [
        "energy, parameter = optimizer.optimize(dimensions=1, function=objective_function)\n",
        "\n",
        "print(f\"\\n=== Adam Optimizer Results ===\")\n",
        "print(f\"Minimized <H> = {energy:.6f}\")\n",
        "print(f\"Optimal parameters: {[round(p, 6) for p in parameter]}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bc4e8d45",
      "metadata": {},
      "source": [
        "### 1.2 SGD (Stochastic Gradient Descent) Optimizer\n",
        "\n",
        "**SGD** is a fundamental optimization algorithm that updates parameters by taking steps proportional to the negative gradient.\n",
        "\n",
        "**Configurable Parameters:**\n",
        "- `step_size` (default: 0.01): Learning rate for parameter updates\n",
        "- `batch_size` (default: 1): Number of samples per batch\n",
        "- `f_tol` (default: 1e-4): Convergence tolerance\n",
        "- `max_iterations`: Maximum number of iterations\n",
        "- `initial_parameters`: Starting parameter values\n",
        "\n",
        "SGD is simpler than Adam and can be effective when you understand your problem well enough to tune the learning rate appropriately."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "491cd22c",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "=== SGD Optimizer Results ===\n",
            "Minimized <H> = -1.748865\n",
            "Optimal parameters: [-5.688733]\n"
          ]
        }
      ],
      "source": [
        "# Configure SGD optimizer\n",
        "sgd_optimizer = cudaq.optimizers.SGD()\n",
        "sgd_optimizer.step_size = 0.05       # Learning rate\n",
        "sgd_optimizer.batch_size = 1         # Stochastic mode\n",
        "sgd_optimizer.max_iterations = 100   # Maximum iterations\n",
        "sgd_optimizer.f_tol = 1e-6           # Convergence tolerance\n",
        "sgd_optimizer.initial_parameters = initial_params\n",
        "\n",
        "# Run optimization\n",
        "sgd_energy, sgd_params = sgd_optimizer.optimize(dimensions=1, function=objective_function)\n",
        "\n",
        "print(f\"\\n=== SGD Optimizer Results ===\")\n",
        "print(f\"Minimized <H> = {sgd_energy:.6f}\")\n",
        "print(f\"Optimal parameters: {[round(p, 6) for p in sgd_params]}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4caa799e",
      "metadata": {},
      "source": [
        "### 1.3 SPSA (Simultaneous Perturbation Stochastic Approximation)\n",
        "\n",
        "**SPSA** is a gradient-free stochastic optimization algorithm that is particularly useful for noisy objective functions (like quantum hardware with shot noise). It approximates gradients using simultaneous perturbations and requires only **2 function evaluations per iteration** regardless of problem dimension.\n",
        "\n",
        "**Configurable Parameters:**\n",
        "- `step_size` (default: 0.3): Evaluation step size for gradient approximation\n",
        "- `gamma` (default: 0.101): Scaling exponent for step size schedule\n",
        "- `max_iterations`: Maximum number of iterations\n",
        "- `initial_parameters`: Starting parameter values\n",
        "\n",
        "**Key Advantage**: SPSA does **not** require gradients, making it ideal for noisy functions and quantum hardware."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "85afa4dd",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "=== SPSA Optimizer Results ===\n",
            "Minimized <H> = -1.748668\n",
            "Optimal parameters: [-5.681724]\n"
          ]
        }
      ],
      "source": [
        "# Configure SPSA optimizer\n",
        "spsa_optimizer = cudaq.optimizers.SPSA()\n",
        "spsa_optimizer.step_size = 0.3       # Evaluation step size\n",
        "spsa_optimizer.gamma = 0.101         # Scaling exponent\n",
        "spsa_optimizer.max_iterations = 100  # Maximum iterations\n",
        "spsa_optimizer.initial_parameters = initial_params\n",
        "\n",
        "# Define gradient-free objective function\n",
        "def spsa_objective(parameter_vector: list[float]) -> float:\n",
        "    \"\"\"\n",
        "    Objective function for gradient-free optimizers like SPSA.\n",
        "    Returns: cost only (no gradient)\n",
        "    \"\"\"\n",
        "    return cudaq.observe(kernel, hamiltonian, parameter_vector).expectation()\n",
        "\n",
        "# Run optimization\n",
        "spsa_energy, spsa_params = spsa_optimizer.optimize(dimensions=1, function=spsa_objective)\n",
        "\n",
        "print(f\"\\n=== SPSA Optimizer Results ===\")\n",
        "print(f\"Minimized <H> = {round(spsa_energy, 6)}\")\n",
        "print(f\"Optimal parameters: {[round(p, 6) for p in spsa_params]}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a36a46fc-b04b-43e3-ac10-3b2ee59c3418",
      "metadata": {},
      "source": [
        "## 2. Third-Party Optimizers\n",
        "\n",
        "CUDA-Q optimizers can work alongside third-party optimization libraries like SciPy. This provides flexibility to use familiar optimization tools while leveraging CUDA-Q's quantum simulation capabilities.\n",
        "\n",
        "The same VQE procedure can be accomplished using SciPy. In this case, a simple cost function is defined and provided as input to the standard SciPy `minimize` function."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "0541fa5c-3c4e-44e7-bd13-52d3b07af50e",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            " message: Optimization terminated successfully.\n",
            " success: True\n",
            "  status: 1\n",
            "     fun: -1.748865011330396\n",
            "       x: [ 5.943e-01]\n",
            "    nfev: 26\n",
            "   maxcv: 0.0\n"
          ]
        }
      ],
      "source": [
        "from scipy.optimize import minimize\n",
        "\n",
        "def cost(theta):\n",
        "\n",
        "    exp_val = cudaq.observe(kernel, hamiltonian, theta).expectation()\n",
        "\n",
        "    return exp_val\n",
        "\n",
        "result = minimize(cost, initial_params ,method='COBYLA', options={'maxiter': 40})\n",
        "print(result)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d73755a8-2c45-41c8-8d9b-d88a0b411869",
      "metadata": {},
      "source": [
        "## 3. Parallel Parameter Shift Gradients"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7b01d93b-0953-4139-8a90-44414327f726",
      "metadata": {},
      "source": [
        "CUDA-Q's `mqpu` backend allows for parallel computation of parameter shift gradients using multiple simulated QPUs. Gradients computed this way can be used in any of the previously discussed optimization procedures.  Below is an example demonstrating how parallel gradient evaluation can be used for a VQE procedure. \n",
        "\n",
        "The parameter shift procedure computes two expectations values for each parameter shifted forwards and backwards. These are used to estimate the gradient contribution for that parameter.\n",
        "\n",
        "The following code defines a function that takes a kernel, a Hamiltonian (spin operator), and the circuit parameters and produces a parameter shift gradient with shift `epsilon`. The first step of the function builds `xplus` and `xminus` , arrays consisting of the shifted parameters. \n",
        "\n",
        "Next, a for loop iterates over all of the parameters and uses the `cudaq.observe_async` to compute the expectation value.  This command also takes `qpu_id` as an in out which specifies the GPU that will be used to simulate the ith QPU.  In the example below, four GPUs (simulated QPUs) are available so the gradient is batched over four devices. \n",
        "\n",
        "The results are saved in the `g_plus` and `g_minus` lists, the elements of which are accessed with commands like  `g_plus[1].expectation()` to compute the finite differences and construct the final gradient. \n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "7508c92a-ff80-4a92-9396-691825b1cabb",
      "metadata": {},
      "outputs": [],
      "source": [
        "import  numpy as np\n",
        "# cudaq.set_target('nvidia', option = 'mqpu')\n",
        "\n",
        "num_qpus = 1\n",
        "epsilon =np.pi/4\n",
        "\n",
        "\n",
        "def batched_gradient_function(kernel, parameters, hamiltonian, epsilon): \n",
        "\n",
        "    # Prepare an array of parameters corresponding to the plus and minus shifts\n",
        "    x = np.tile(parameters, (len(parameters),1))\n",
        "    xplus = x + (np.eye(x.shape[0]) * epsilon)\n",
        "    xminus = x - (np.eye(x.shape[0]) * epsilon)\n",
        "\n",
        "    g_plus = []\n",
        "    g_minus = []\n",
        "    gradient = []\n",
        "\n",
        "    qpu_counter = 0 # Iterate over the number of GPU resources available\n",
        "    \n",
        "    \n",
        "    for i in range(x.shape[0]): \n",
        "\n",
        "        g_plus.append(cudaq.observe_async(kernel,hamiltonian, xplus[i], qpu_id = qpu_counter%num_qpus))\n",
        "        qpu_counter += 1 \n",
        "\n",
        "        g_minus.append(cudaq.observe_async(kernel, hamiltonian, xminus[i], qpu_id = qpu_counter%num_qpus))\n",
        "        qpu_counter += 1 \n",
        "        \n",
        "    # Use the expectation values to compute the gradient    \n",
        "    gradient = [(g_plus[i].get().expectation() - g_minus[i].get().expectation()) / (2*epsilon) for i in range(len(g_minus))]\n",
        "\n",
        "    return gradient\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9a317073-0e29-40d1-bf09-d168e7227f14",
      "metadata": {},
      "source": [
        "This function can be used in a VQE procedure as presented below. The `batched_gradient_function` is used to evaluate the gradient at each optimization step. This objective function returns the cost and gradient at the current parameter values and can be used with any SciPy optimizer that uses gradients (like L-BFGS-B)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "6313f4ed-a620-4982-b587-42e19267ef84",
      "metadata": {},
      "outputs": [],
      "source": [
        "def objective_function(parameter_vector,\n",
        "                       hamiltonian=hamiltonian,\n",
        "                       kernel=kernel,\n",
        "                       epsilon=epsilon):\n",
        "    \"\"\"\n",
        "    Objective function for VQE with parallel parameter shift gradients.\n",
        "    Computes both cost and gradient at the current parameter values.\n",
        "    \"\"\"\n",
        "    # Compute cost at current parameters\n",
        "    cost = cudaq.observe(kernel, hamiltonian, parameter_vector).expectation()\n",
        "    \n",
        "    # Compute gradient at current parameters using parallel parameter shift\n",
        "    gradient_vector = batched_gradient_function(kernel, parameter_vector, hamiltonian, epsilon)\n",
        "\n",
        "    return cost, gradient_vector\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "edd702a6-213d-469e-861b-9beee6207284",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Run VQE optimization with parallel parameter shift gradients\n",
        "result_vqe = minimize(objective_function, initial_params, method='L-BFGS-B', jac=True, tol=1e-8, options={'maxiter': 50})\n",
        "\n",
        "print(\"\\n=== VQE with Parallel Parameter Shift Gradients ===\")\n",
        "print(f\"Optimized energy: {result_vqe.fun:.6f}\")\n",
        "print(f\"Optimal parameters: {result_vqe.x}\")\n",
        "print(f\"Number of iterations: {result_vqe.nit}\")\n",
        "print(f\"Success: {result_vqe.success}\")"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}
