Hello everyone! I am researching the use of different numerical methods for calculating A−1 when with A a square matrix of a certain size (eg. 3x3). I would like to share my methods with you and seek your opinions.

Method 1: Use of the formula A−1=Com(A)T/det(A)

Method 2: Resolution by the Gauss algorithm of the unknown system X∈M10,1​(R), where (a1​,...,a10​) are arbitrary (S)A⋅X=​a1​...a10​​​

Method 3: Resolution of system (S) by triangularizing A.

Method 4: Resolution of system (S) using the QR form from the course on Euclidean spaces (we saw and programmed it in td-tp).

We don't use ready-made system resolution functions in Python, reduction function, the determinant calculation function.

How to compare different methods: we compare the number of basic operations (+,-,x,/) of each method. As we don't know how many calculations inside the det(), eigenvects(), etc, we don't allow using them and only use the most basic calculation.

Hi, I'm trying to solve the KdV equation with the Crank-Nikolson Scheme and I'm trying to follow the method in this document (pg4). I am getting confused on how to iterate my loops because of all of the different indices and keep track of values etc. If anyone could give any advice, that'd be wonderful. Thank you! :)

Hello guys,

As a first disclaimer, and maybe as a first question, I know very little about optimization, so probably the doubts in this post will come from my lack of knowledge. Is there any standard reference to introduce me in the different optimization methods?

The main point of this post is that I'm performing a minimization of an objective function using the scipy.optimize.minimize package.If I use BFGS it says it performs 0 iterations, if I use Nelder-Mead, it iterates but the final value is the initial condition. I can't figure if it is code problem or concepts problem. Here is my code if it helps:

import numpy as np
import scipy.sparse as sp
import scipy.optimize as so
from import_h5py import K_IIDF, K_IBD, K_IBD_T, K_BBD
from grids_temperaturas import C_RD_filtrada
from leer_conductors import conductividades

def construct_KR(conduct_dict):
    """
    Convierte un diccionario de conductividades a matriz sparse.
    
    Args:
        conduct_dict: Diccionario donde las claves son tuplas (i,j) y los valores son conductividades.
        
    Returns:
        sp.csc_matrix: Matriz sparse en formato CSC.
    """
    if not conduct_dict:
        return sp.csc_matrix((0, 0))
    
    # Encontrar dimensiones máximas
    max_i = max(k[0] for k in conduct_dict.keys())
    max_j = max(k[1] for k in conduct_dict.keys())
    n = max(max_i, max_j)  # Asumimos matriz cuadrada
    
    # Preparar datos para matriz COO
    rows, cols, data = [], [], []
    for (i, j), val in conduct_dict.items():
        rows.append(i-1)  # Convertir a 0-based index
        cols.append(j-1)
        data.append(val)
    
    # Crear matriz simétrica (sumando transpuesta)
    K = sp.coo_matrix((data + data, (rows + cols, cols + rows)), shape=(n, n))
    row_sums = np.array(K.sum(axis=1)).flatten()
    Kii = sp.diags(row_sums, format='csc')
    K_R = Kii-K

    boundary_nodes = []
    
    with open('bloque_nastran3_acase.BOUNDARY_CONDS.data', 'r') as f:
        for line in f:
            # Busca líneas que definen temperaturas de contorno
            if line.strip().startswith('T') and '=' in line:
                # Extrae el número de nodo (ej. 'T37' -> 37)
                node_str = line.split('=')[0].strip()[1:]  # Elimina 'T' y espacios
                try:
                    node_num = int(node_str)
                    boundary_nodes.append(node_num - 1)  # Convertir a 0-based
                except ValueError:
                    continue

        """
    Reordena la matriz y vector según nodos libres y con condiciones de contorno.
    
    Args:
        sparse_matrix: Matriz de conductividad (Kii - K) en formato sparse
        temperature_vector: Vector de temperaturas
        boundary_file: Ruta al archivo .BOUNDARY_CONDS
        
    Returns:
        tuple: (matriz reordenada, vector reordenado, número de nodos libres)
    """
    # Leer nodos con condiciones de contorno del archivo
    constrained_nodes = boundary_nodes
    size = K_R.shape[0]
    
    # Verificar que los nodos de contorno son válidos
    for node in constrained_nodes:
        if node < 0 or node >= size:
            raise ValueError(f"Nodo de contorno {node+1} está fuera de rango")

    # Crear máscaras booleanas
    bound_mask = np.zeros(size, dtype=bool)
    bound_mask[constrained_nodes] = True
    free_mask = ~bound_mask

    # Índices de nodos libres y con condición de contorno
    free_nodes = np.where(free_mask)[0]
    constrained_nodes = np.where(bound_mask)[0]

    # Nuevo orden: primero libres, luego con condición de contorno
    new_order = np.concatenate((free_nodes, constrained_nodes))

    num_free = len(free_nodes)
    num_constrained = len(constrained_nodes)
    # Reordenar matriz y vector
    K_ordered = sp.csc_matrix(K_R[new_order][:, new_order])

    K_IIRF = K_ordered[:num_free, :num_free]

    K_IBR = K_ordered[:num_free,:num_constrained]

    K_IBR_T = K_IBR.transpose()

    K_BBR = K_ordered[:num_constrained,:num_constrained]

    return K_IIRF,K_IBR,K_IBR_T,K_BBR

#K_IIRF_test,_,_,_ = construct_KR(conductividades)
#resta = K_IIRF - K_IIRF_test
#print("Norma de diferencia:", np.max(resta))

## Precalculos

K_IIDF_K_IBD = sp.linalg.spsolve(K_IIDF, K_IBD)
K_IIDF_KIBD_C_RD = C_RD_filtrada @ sp.linalg.spsolve(K_IIDF, K_IBD)

def calcular_epsilon_CMF(cond_vector):
    """
    Calcula epsilon_CMF según la fórmula proporcionada.

    Args:
        epsilon_MO: Matriz de errores MO de tamaño (n, n_b).
        epsilon_MT: Matriz de errores MT de tamaño (n_r, n_b).
        
    Returns:
        epsilon_CMF: Valor escalar resultante.
    """
    nuevas_conductividades = cond_vector.tolist()
    nuevo_GL = dict(zip(vector_coordenadas, nuevas_conductividades))

    K_IIRF,K_IBR,K_IBR_T,K_BBR = construct_KR(nuevo_GL)


    #epsilon_MQ = sp.linalg.spsolve(K_BBD,sp.csc_matrix(-K_IBD_T @ sp.linalg.spsolve(K_IIDF, K_IBD) + K_BBD) - (-K_IBR_T @ sp.linalg.spsolve(K_IIRF, K_IBR) + K_BBR )) 

    epsilon_MQ = sp.linalg.spsolve(K_BBD,((-K_IBD_T @ K_IIDF_K_IBD + K_BBD) - (-K_IBR_T @ sp.linalg.spsolve(K_IIRF, K_IBR) + K_BBR )))


    epsilon_MT = sp.linalg.spsolve(K_IIRF,K_IBR) - K_IIDF_KIBD_C_RD
    # Suma de cuadrados (usando .power(2) y .sum() para sparse)
    sum_MQ = epsilon_MQ.power(2).sum()
    sum_MT = epsilon_MT.power(2).sum()

    
    # Raíz cuadrada del total
    epsilon_CMF = np.sqrt(sum_MQ + sum_MT)
    
    return epsilon_CMF

def debug_callback(xk):
    print("Iteración, error:", calcular_epsilon_CMF(xk))


cond_opt = so.minimize(
    calcular_epsilon_CMF,
    vector_conductividades,
    method='Nelder-Mead',
    options={'maxiter': 100, 'xatol': 1e-8},
    callback=debug_callback
)
print("epsilon_CMF:", cond_opt)

The idea is that, with each iteration, the cond_vector changes and then it generates new matrices with the construct_KR, so it recalculates epsilon_CMF

We are currently learning Matrix equations in pre cal and I have come to realize that it takes a good 5 minutes to solve one equation. There has to got be a quicker and more efficient way of solving these problems.

I’d like to know the historical process behind two mathematical/numerical methods:

• Polar decomposition (factorizing a matrix).
• Newton–Schulz iteration.

My question isn’t just who first wrote them down, but how they were invented:

• What problems or contexts led Autonne (or others) to polar decomposition? Was it geometric (analogy to complex polar form), mechanical (deformation gradient = rotation × stretch), or theoretical?
• How did Schulz’s idea emerge? Was it a response to early computational limitations, or a mathematical curiosity later applied to matrices?

I’d love to understand what kind of analogies, problems, or constraints guided these mathematicians — essentially, how they thought their way into discovering these methods, not just the final result. I’d appreciate a timeline, the key figures/papers, and especially what the inventors were trying to achieve at the time.

Hey! How should I solve h(t) from this system of these equations (if it's even possible, sorry if it's a dumb question) https://i.gyazo.com/f4df0be6db25e59cb8b10c3610caf3ba.png numerically(euler's method or something like that...?)

Initial conditions being something like this : h(0) = 0, h'(0) = 0 h''(0) = 0

Hi, i'm relatively new to both python and math (I majored in history something like a year ago) so i get if the problem i'm about to ask help for sounds very trivial.

My code has started running super slow out of nowhere, i was literally running it in 30 seconds, despite the multiple nested loops that calculated 56 million combinations, it was relatively ok even with a very computationally heavy grid search for my parameters. I swear, i went to get coffee, did not even turn down the pc, from one iteration to the other now 30 minutes of waiting time. Mind you, i have not changed a single thing

(these are three separate pi files, just to illustrate the process I'm going through)

FIRST FILE:

std = np.linalg.cholesky(matrix)

part = df['.ARTKONE returns'] + 1

ψ = np.sqrt(np.exp(np.var(part) - 1))
emp_kurtosis = 16*ψ**2 + 15*ψ**4 + 6*ψ**6 + ψ**8
emp_skew = 3*ψ + ψ**3

intensity = []
jump_std = []
brownian_std = []

for λ in np.linspace(0,1,100): 
    for v in np.linspace(0,1,100):
        for β in np.linspace(0,1,100):
            ξ = np.sqrt(np.exp(λ*v**2 + λ*β**2) - 1)
            jump_kurtosis = 16*ξ**2 + 15*ξ**4 + 6*ξ**6 + ξ**8     
            jump_skew = 3*ξ + ξ**3
            if np.isclose(jump_kurtosis,emp_kurtosis, 0.00001) == True and np.isclose(emp_skew,jump_skew, 0.00001) == True:
                print(f'match found for: - intensity: {λ} -- jump std: {β} -- brownian std: {v}') 

SECOND FILE:

df_3 = pd.read_excel('paraameters_values.xlsx')
df_3.drop(axis=1, columns= 'Unnamed: 0', inplace=True)

part = df['.ARTKONE returns'] + 1

mean = np.mean(part)
ψ = np.sqrt(np.exp(np.var(part) - 1))
var_psi = mean * ψ

for i in range(14):

    λ = df_3.iloc[i,0]
    β = df_3.iloc[i,1]
    v = df_3.iloc[i,2]

    for α in np.linspace(-1,1,2000):
        for δ in np.linspace(-1,1,2000):
            exp_jd_r = np.exp(δ +λ - λ*(np.exp(α - 0.5 * β **2)) + λ*α + λ*(0.5 * β **2))
            var_jd_p =  (np.sqrt(np.exp(λ*v**2 + λ*β**2) - 1)) * exp_jd_r **2 
            if np.isclose(var_jd_p, var_psi, 0.0001) == True and np.isclose(exp_jd_r, mean, 0.0001) == True:
                print(f'match found for: - intensity: {λ} -- jump std: {β} -- brownian std: {v} -- delta: {δ} -- alpha: {α}')

FUNCTIONS: because (where psi is usally risk tolerance = 1, just there in case i wanted a risk neutral measure)

def jump_diffusion_stock_path(S0, T, μ, σ, α, β, λ, φ):
    n_j = np.random.poisson(λ * T)
    μj = μ - (np.exp(α + 0.5*β**2) -1) * λ *φ + ((n_j * np.log(np.exp(α + 0.5*β**2)))/T)
    σj = σ**2 + (n_j * β **2)/T 
    St = S0 * np.exp(μj * T - σj * T * 0.5 + np.sqrt(σj * T) * np.random.randn())
    return St
def geometric_brownian_stock_path(S0, T, μ, σ):
    
    St = S0 * np.exp((μ-(σ**2)/2)*T + σ * np.sqrt(T) * np.random.randn())
    return St

I know this code looks ghastly, but given it was being handled just fine, and all of a sudden it didn't, i cannot really explain this. I restarted the pc, I checked memory and cpu usage (30, and 10% respectively) using mainly just two cores, nothing works.
i really cannot understand why, it is hindering the progression of my work a lot because i rely on being able to make changes quickly as soon as i see something wrong, but now i have two wait 30 minutes before even knowing what is wrong. One possible issue is that these files are in folders where multiple py files call for the same datasets, but they are inactive so this should not be a problem.

:there's no need to read this second part, but i put it in if you're interested

THE MATH: I'm trying to define a distribution for a stochastic process in such a way that it resembles the empirical distribution observed in the past for this process (yes the data i have is stationary), to do this i'm trying to build a jump diffusion process (lognormal, poisson, normally distributed jump sizes). In order for this jump diffusion process to match my empirical distribution i created two systems of equations: one where i equated the expected value of the standard brownian motion with the one of the jump diffusion, and did the same for the expected values of their second moments, and a second where i equated the kurtosis of the empirical distribution to the standardised fourth moment of the jump diffusion, and the skew of the empirical to the third standardised moment of the jump diffusion.
Since i am too lazy to go and open up a book and do it the right way or to learn how to set up a maximum likelihood estimation i opted for a brute gride search.
Why all this??
i'm working on inserting alternative assets in an investment portfolio, namely art, in order to do so with more advance techniques, such as CVaR or the jacobi bellman dynamic programming approach, i need to define the distribution of my returns, and art returns are very skewed and and have a lot of kurtosis, simply defining their behaviour as a lognormal brownian motion with N(mean, std) would cancel out any asymmetry which characterises the asset.

thank you so much for your help, hope you all have a lovely rest of the day!