This became a little bit convoluted but it seems to be working.

For each pair of points azimuth is calculated then translated into directions (N, E ,S or W)

Change layer and field names to match your data:

dm = QgsProject.instance().mapLayersByName('Distance matrix')[0]
pointlyr1 = QgsProject.instance().mapLayersByName('ns_point')[0] #Input point layer in distance matrix
pointlyr1_idfield = 'ogc_fid'
pointlyr2 = QgsProject.instance().mapLayersByName('js_point')[0] #Target 
pointlyr2_idfield = 'ogc_fid'

new_field_name = 'Direction' #Add this text field to distance matrix before executing the code

dms = [[int(f['InputID']),int(f['TargetID'])] for f in dm.getFeatures()] #List all input and target ids
p1dict = {int(f[pointlyr1_idfield]):f.geometry() for f in pointlyr1.getFeatures()} #Create a dictionary of id and geometry for input
p2dict = {int(f[pointlyr2_idfield]):f.geometry() for f in pointlyr2.getFeatures()} #Same for target

directions = [] #List to hold directions (N,E,S,W)

for id1, id2 in dms:
    #Function to calculate azimuth and then translate this to direction
    az = p1dict[id1].asPoint().azimuth(p2dict[id2].asPoint())
    def calcdir(a):
        if -45<a<45:
            v = 'N'
        elif 45<=a<135:
            v = 'E'
        elif 135<=a<=180 or -180<=a<=-135:
            v = 'S'
        else:
            v = 'W'
        return v
    directions.append(calcdir(az))
#Write directions list to a column in Distance matrix layer
with edit(dm):
    for feat, direct in zip(dm.getFeatures(), directions):
        #print(feat.id(), direct)
        feat.setAttribute(feat.fieldNameIndex(new_field_name), direct) 
        dm.updateFeature(feat)

This became a little bit convoluted but it seems to be working.

For each pair of points azimuth is calculated then translated into directions (N, E ,S or W)

Change layer and field names to match your data:

dm = QgsProject.instance().mapLayersByName('Distance matrix')[0]
pointlyr1 = QgsProject.instance().mapLayersByName('ns_point')[0] #Input point layer in distance matrix
pointlyr1_idfield = 'ogc_fid'
pointlyr2 = QgsProject.instance().mapLayersByName('js_point')[0] #Target 
pointlyr2_idfield = 'ogc_fid'

new_field_name = 'Direction' #Add this text field to distance matrix before executing the code

dms = [[int(f['InputID']),int(f['TargetID'])] for f in dm.getFeatures()] #List all input and target ids
p1dict = {int(f[pointlyr1_idfield]):f.geometry() for f in pointlyr1.getFeatures()} #Create a dictionary of id and geometry for input
p2dict = {int(f[pointlyr2_idfield]):f.geometry() for f in pointlyr2.getFeatures()} #Same for target

directions = [] #List to hold directions (N,E,S,W)

for id1, id2 in dms:
    #For each pair of points calculate azimuth (-180 - +180)
    az = p1dict[id1].asPoint().azimuth(p2dict[id2].asPoint())
    def calcdir(a):
        #Translate azimuth to compass direction
        if -45<a<45:
            v = 'N'
        elif 45<=a<135:
            v = 'E'
        elif 135<=a<=180 or -180<=a<=-135:
            v = 'S'
        else:
            v = 'W'
        return v
    directions.append(calcdir(az))
#Write directions list to a column in Distance matrix layer
with edit(dm):
    for feat, direct in zip(dm.getFeatures(), directions):
        #print(feat.id(), direct)
        feat.setAttribute(feat.fieldNameIndex(new_field_name), direct) 
        dm.updateFeature(feat)